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Abstract 


Recent  advances  in  computational  plasma  physics  have  been 
made  by  using  a  fluid  particle  approach  to  simulate  the  time 
dependent,  compressible  MHD  equations.  Inviscid  fluids  in 
periodic  geometries  have  been  successfully  modelled  by  other 
researchers  [2], [7]  with  such  a  particle  method.  The 
development  of  a  method  for  the  computer  modelling  of 
viscous,  compressible  fluids  in  nonperiodic  geometries  using 
a  fluid  particle  approach  is  described  in  the  following 
pages . 

The  fluid  equations  are  first  written  in  a  Lagrangian 
form,  which  is  conservative  in  a  Lagrangian  sense.  The 
fluid  particle  approach  then  represents  the  fluid  by  a  large 
number  of  "particles",  which  are  not  true  point  particles 
but  rather  have  a  finite  spatial  extent.  Each  finite 
particle  has  a  velocity  and  temperature,  so  that  the  fluid 
quantities  of  velocity  and  temperature  are  obtained  by 
accumulating  the  particle  quantities  on  a  fixed 
computational  grid.  A  small  scale  smoothing  is  applied  to 
the  fluid  velocity  and  temperature  to  reduce  noise  levels 
caused  by  random  particle  motion.  The  finite  particles  are 
given  a  Gaussian  shaped  mass  density,  so  that  the  fluid  mass 
density  is  obtained  by  summing  the  mass  densities  of  nearby 
particles  using  a  nearest  grid  point  (NGP)  approach.  The 
fluid  equations  are  then  finite  differenced  on  the  grid,  and 
the  Lagrangian  derivatives  appearing  in  the  fluid  equations 
are  then  used  to  advance  the  particle  velocities  and 


v 


temperatures.  Particle  positions  are  changed  according  to 
particle  velocities.  A  new  method  of  assigning  the  pressure 
gradient  force  to  the  particles  is  developed,  which  resulted 
in  a  very  significant  reduction  in  the  undesirable  random 
particle  motion  that  occurred  with  the  inviscid  simulations 
done  by  other  researchers  [2],  [7].  The  flow  of  a  turbulent 
viscous  fluid  through  a  two  dimensional  channel  is  described 
analytically  and  experimentally,  and  a  fluid  particle  method 
for  simulating  the  channel  is  then  developed.  Several 
special  considerations,  required  at  the  open  flow 
boundaries,  are  described.  The  code  results  are  compared  to 
the  analytic  and  experimental  results.  It  is  found  that  the 
multistreaming  terms  (which  prevent  particles  from  having 
velocities  and  temperatures  much  different  from  the  fluid 
velocity  and  temperature,  as  used  in  [2])  must  be  used 
carefully  in  the  case  of  channel  flow.  The  source  of  the 
random  particle  noise  is  also  discussed. 

A  magnetically  stabilized  gas  discharge,  which  was 
being  experimentally  tested  at  the  U  of  A  [4],  and  was 
modelled  numerically  by  Antoniuk  in  [4], [16]  using  a  conven¬ 
tional  fluid  approach,  is  also  modelled  and  the  results 
compared  to  Antoniuk’ s  experimental  and  computational 
results . 
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CHAPTER  I 


Introduction 

Much  of  the  nature  of  the  physical  world  can  be  described  by 
partial  differential  equations.  Before  the  advent  of  the 
digital  computer  only  a  limited  number  of  these  equations 
had  been  solved  to  the  extent  of  yielding  useful  physical 
solutions.  However,  in  the  past  twenty  years  many  solutions 
of  these  equations  have  been  obtained  by  computer 
simulation . 

In  modelling  any  set  of  partial  differential  equations 
on  a  computer  it  is  necessary  to  break  the  continuous  nature 
of  the  partial  differential  equations  into  some  discrete 
form.  One  way  of  doing  this  is  the  finite  difference 
method.  In  the  case  where  the  differential  equations 
involve  the  spatial  and  temporal  variation  of  several  quan¬ 
tities,  then  finite  differencing  begins  by  considering  these 
quantities  at  only  a  finite  number  of  points  in  time  and 
space.  These  points  form  a  temporal,  spatial  mesh  or  grid. 
The  derivatives  appearing  in  the  partial  differential 
equations  are  then  approximated  on  the  grid  by  the  differ¬ 
ence  between  the  quantities  at  neighboring  grid  points. 

Using  these  approximated  derivatives  in  the  equations,  the 
values  of  the  variables  involved  in  the  equations  can  be 
evaluated  at  each  grid  point. 

In  the  area  of  fluid  dynamics  there  have  been  two  basic 
approaches  to  applying  the  finite  difference  method  to  the 
governing  partial  differential  equations  of  fluid  motion. 
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The  Eulerian  codes  (code  is  a  term  meaning  computer  program) 
use  a  f ixed  spatial  finite  difference  grid,  across  which  the 
fluid  streams  by.  The  grid  is  called  Eulerian  because  it  is 
stationary,  so  that  the  derivatives  evaluated  directly  on 
the  grid  are  partial  or  Eulerian  derivatives.  The  other 
scheme  involves  moving  the  spatial  grid  with  the  fluid, 
resulting  in  a  Lagrangian  code.  The  term  Lagrangian  refers 
to  the  fact  that  the  spatial  grid  distorts  with  the  fluid, 
and  the  derivatives  on  the  grid  are  Lagrangian  derivatives 
(also  known  as  convective,  substantive  or  total 
derivatives) . 

One  of  the  most  severe  drawbacks  of  the  Eulerian  method 
is  the  difficulty  found  in  dealing  correctly  with  the 
advective  terms  (L i { v , 9A/dx j } )  in  the  total  derivatives 
dA/dt  appearing  in  the  fluid  equations.  Instabilities 
producing,  for  example,  negative  densities  can  arise.  These 
can  be  dealt  with  by  such  elaborate  methods  as  flux 
corrected  transport  (FCT)  [15],  but  this  adds  complexity  as 
well  as  expense.  As  an  alternative,  the  Lagrangian  approach 
eliminates  the  advective  derivative  problems;  however,  in 
two  and  three  dimensions,  the  spatial  finite  difference  grid 
becomes  so  entangled  that  the  method  breaks  down.  As  a 
result,  in  more  than  one  dimension  the  Eulerian  method  is 
generally  used. 

In  plasma  physics,  the  laws  governing  a  plasma  can  be 
written  in  a  form  similar  to  the  normal  fluid  dynamics 
equations.  For  a  plasma  which  is  assumed  to  be 
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quasi-neutral  (i.e.  charge  density  p= 0),  the 
magnetohydrodynamic  (MHD)  equations  result.  These  are 
essentially  an  extension  of  the  normal  fluid  equations  to  a 
conducting  fluid,  with  the  addition  of  Maxwell’s 
electromagnetic  equations.  Alternatively,  the  plasma  may  be 
treated  as  a  multi-species  fluid;  in  which  case,  the 
electrons,  ions,  and  neutrals  are  treated  as  separate, 
interacting  fluids,  each  obeying  equations  with  a  similar 
form  to  the  MHD  equations.  The  Eulerian  or  Lagrangian 
finite  difference  method  can  be  applied  to  the  sets  of 
equations  resulting  from  the  MHD  or  multi-species 
treatments;  however,  in  addition  to  the  finite  difference 
difficulties  mentioned  above,  these  methods  yield  solutions 
to  a  plasma  which  has  been  assumed  to  obey  a  fluid  picture. 
Many  real  plasmas,  and  most  nonlinear  plasma  phenomena,  do 
not  obey  the  fluid  type  plasma  equations  and  so  the  Eulerian 
or  Lagrangian  fluid  finite  difference  schemes  cannot  be  used 
with  accuracy. 

To  overcome  this  difficulty  of  plasma  simulations,  one 
might  consider  dealing  with  individual  electrons  and  ions, 
since  their  equations  of  motion  are  quite  simple  and  one 
could  solve  self-consistently  for  the  electric  and  magnetic 
fields  based  on  the  particle  positions  and  velocities. 
Unfortunately,  the  number  of  particles  involved  in  any 
practical  system  is  far  beyond  the  computational  and  memory 
capabilities  of  even  the  most  advanced  supercomputers.  In 
the  1960's  though,  a  simulation  technique  involving  what  has 


been  called  "finite  particles"  arose  [5], [6].  In  most 
plasma  simulations,  these  finite  particles  obey  the  same 
equations  as  a  single  real  particle.  In  addition  to  having 
a  position  and  velocity  these  macroparticles  may  also  carry 
a  charge,  but  since  they  are  finite,  this  charge  is  spread 
out  over  a  finite  volume  in  space.  Since  each  macroparticle 
represents  a  large  number  of  real  particles,  the  number  of 
macroparticles  required  to  model  a  realistic  geometry  is 
within  the  limits  of  today's  computers. 

Earlier  finite  particle  plasma  simulations  advanced  the 
finite  particles  explicitly  in  time.  An  explicit  time  dif¬ 
ference  scheme  evaluates  the  values  of  the  variables  at  time 
(n+1)At  based  solely  on  their  values  at  previous  time  steps. 
Explicit  methods  are  conditionally  stable:  they  must  satisfy 
the  condition  that  the  time  step  At  must  be  smaller  than  the 
time  it  takes  for  the  fastest  phenomena  to  cross  one  grid 
spacing.  This  requirement  is  called  the  Courant  or 
Courant-Fr iedr ichs-Lewy  (CFL)  stability  criterion.  For 
plasma  simulations  this  is  an  extremely  severe  limitation 
since  the  fastest  time  scale  is  determined  by  either  the  ion 
plasma,  electron  cyclotron,  or  photon  frequency,  all  of 
which  usually  yield  time  scales  which  are  many  orders  of 
magnitude  below  the  time  scale  of  interest.  This  places  the 
explicit  simulation  of  most  plasmas  beyond  the  computation 
times  available. 

In  recent  years,  the  Courant  limit  has  been  removed  by 
the  development  of  implicit  time  advancement  methods, 
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implicit  meaning  that  the  method  of  advancing  the  quantities 
at  time  nAt  is  dependent  on  the  values  one  step  ahead  in 
time  at  (n+1)At.  Clearly  this  requires  some  sophisticated 
computation,  and  indeed  implicit  particle  simulations  are 
presently  largely  limited  to  one  and  two  dimensions  because 
of  this.  Some  recent  hybrid  codes  have  combined  a  finite 
particle  method  (to  deal  with  the  nonfluid  electrons  and 
ions)  with  a  multifluid  FCT  Eulerian  scheme  (to  take  care  of 
the  cold  background  electrons  and  ions)  [10]. 

The  particle  plasma  methods  in  general  are  neither 
cheap  nor  simple,  and  for  plasmas  which  can  be  assumed  to 
obey  the  fluid  type  equations,  or  which  are  three  dimen¬ 
sional,  it  is  much  easier  and  cheaper  to  use  an  Eulerian 
finite  difference  scheme.  However,  to  remove  the  advective 
derivative  difficulty  associated  with  these  codes  (mentioned 
earlier),  it  is  possible  to  employ  a  modified  finite 
particle  approach  to  the  fluid  equations,  which  is  cost 
competitive  with  the  Eulerian  codes.  The  application  of 
such  a  fluid  particle  method  to  viscous  fluids  was  the 
purpose  of  this  work.  The  fluid  particle  method  applies  a 
particle  approach  to  the  fluid  equations,  rather  than  to  the 
particle  equations  which  the  particle  plasma  methods  use. 

When  applying  a  finite  particle  method  to  the  fluid 
equations,  it  is  necessary  to  write  the  equations  in  their 
Lagrangian  form.  In  this  way  the  equations  are  written  as 
dA/dt=...  so  that  the  quantity  A  can  be  carried  by  a  finite 
particle.  The  change  (dA/dt)At  is  then  very  naturally 
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associated  with  the  change  in  A  as  the  finite  particle  moves 
from  position  x  to  position  x  =  x+Y.At .  For  the  quasi-neutral 
or  neutral  fluid  equations,  rather  than  giving  a 
macroparticle  a  spatial  charge  density  it  is  more  convenient 
for  them  to  have  a  spatial  mass  density.  The  fluid  is  thus 
represented  by  a  large  number  of  macroparticles  which  may  be 
thought  of  as  having  a  cloud  of  mass  or  a  mass  density 
profile  with  origin  at  the  macroparticle's  center.  In 
addition  to  having  a  position,  velocity,  and  mass  density, 
the  particles  may  carry  a  temperature  or  any  other  quantity 
which  has  as  its  governing  equation:  dA/dt=...  .  This 
approach  to  the  fluid  equations  removes  the  serious 
advective  difficulties  found  in  traditional  Eulerian  finite 
difference  codes  by  having  the  macroparticles  carry  the 
Lagrangian  derivatives.  In  order  to  maintain  a  spatial  res¬ 
olution  of  approximately  one  grid  spacing,  the 
macroparticles  are  given  a  size  of  order  the  same  size  as  a 
grid  cell.  Two  common  finite  particle  shapes  are  shown  in 
figure  1.1.  From  the  locations  of  the  macroparticles  the 
density  of  the  fluid  can  be  evaluated,  and  all  fluid  quan¬ 
tities  can  be  evaluated  at  the  Eulerian  grid  points  as  just 
an  average  of  the  temperatures  and  velocities  of  the 
macroparticles  inside  a  cell.  Electric  and  magnetic  field 
quantities  and  current  densities  are  not  carried  by  the 
particles,  but  can  be  calculated  on  the  grid  by  conventional 
finite  difference  methods  [2], [7]. 
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Figure  1.1  Square  and  guassian  particle  mass  density  shapes 
are  two  density  profiles  often  assigned  to 
particles 
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Previous  finite  particle  approaches  to  the  fluid 
equations  have  considered  nonviscous  (ideal)  fluids  in 
periodic  geometries  [2],  [7],  (The  term  particle  will  from 
now  on  be  used  to  mean  macroparticle  [finite  particle]  as  is 
conventional  in  the  literature.)  Lebouef,  Tajima,  and 
Dawson  in  [2]  use  a  nearest  grid  point  (NGP)  method  to 
assign  the  particles  to  a  cell.  With  the  nearest  grid  point 
method  a  particle's  position  is  approximated  as  being  at  the 
middle  of  the  cell  that  the  particle's  center  is  in  when 
evaluating  fluid  quantities.  In  addition,  all  particles 
whose  centers  are  in  this  cell  are  pushed  with  the  same 
force,  which  is  calculated  using  conventional  finite 
differencing  of  the  fluid  quantities.  The  NGP  approach  is 
computationally  quite  simple,  but  it  leads  to  excessive 
random  particle  motion  (noise)  in  the  results.  As  an 
alternative,  area  weighted  or  particle  in  cell  (PIC)  methods 
[5],  [7],  which  use  square  or  rectangular  particles,  assign 
particle  quantities  to  grids  by  considering  the  particle's 
area  overlap  in  each  cell.  Forces  are  assigned  to  particles 
in  a  similar  way.  The  PIC  schemes  reduce  noise  levels,  and 
with  second  order  PIC  methods  [8]  noise  levels  and  the 
numerical  diffusion  of  quantities  can  be  reduced  further. 
However,  considerable  expense  and  complexity  are  added  by 
employing  area  weighting  methods.  Subtracted  dipole  methods 
(SUD)  [7]  are  also  considerably  more  complex  than  the  NGP 


method . 
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The  intent  of  this  project  has  been  to  develop  a  NGP 
particle  method  for  simulating  compressible,  viscous  fluids 
using  the  full  instantaneous  fluid  (MHD)  equations.  A 
modified  NGP  approach  was  developed  which  yielded  levels  of 
random  particle  motion  which  were  much  lower  than  in 
previous  NGP  methods.  In  adding  viscosity  to  the  fluid 
equations,  the  assumption  of  a  spatially  periodic  system  was 
no  longer  possible  and  special  considerations  had  to  be 
given  to  open  flow  boundaries.  Naturally  the  fast  Fourier 
transform  techniques  used  in  evaluating  the  Vp  force  for  a 
periodic  system  [7]  cannot  be  used  with  accuracy  because  of 
leakage  and  aliasing  errors  introduced  by  the  nonperiodicity 
involved  [22,  p.  278],  Instead,  pressure  was  obtained  from 
the  ideal  gas  law  in  which  the  fluid  density  was  evaluted  as 
an  explicit  sum  of  the  macroparticle  density  contributions. 

A  modified  NGP  approach  was  then  used  to  assign  the  Vp  force 
to  the  particles.  At  open  flow  boundaries,  special  tech¬ 
niques  were  necessary  to  maintain  a  steady  flow  rate. 

By  removing  the  periodic  boundary  conditions  and  using 
the  fully  viscous  fluid  equations,  many  real  problems  can  be 
dealt  with  using  a  particle  approach.  In  addition,  con¬ 
siderable  simplicity  is  gained  over  other  methods  by  using 
particles  to  deal  with  the  total  derivatives. 


CHAPTER  II 


Mathematical  and  Computational  Model 

2.1  Mathematical  model 

In  order  to  use  a  particle  approach  to  model  the  fluid 
equations,  the  partial  differential  equations  which  govern  a 
Newtonian  fluid  must  be  written  in  a  Lagrangian  form.  In 
first  defining  the  form  of  these  equations  it  is  most  con¬ 
venient  to  consider  a  neutral  fluid.  The  equations  for  a 
MHD  fluid  (see  Chapter  5)  or  a  multifluid  plasma  are  very 
similar  and  are  easily  formed  from  the  neutral  fluid 
equations . 

Defining  the  velocity  of  the  fluid  as  v(jr,t)  and  its 
density  by  p(_r,t),  then  the  momentum  equation  can  be  written 
as  [ 1 ,  p.  51 ] 

pdv/dt  =  B  +  V  •  g  (2.1) 

where  B  contains  all  body  forces  acting  on  the  fluid,  g  is 
the  stress  tensor: 

g  =  [-p+ ( /3-2p/3  )V*  v  ]I  +  2pS 

and  the  total  derivative  is  dv/dt  =  9v/3t  +  v*Vv.  Here,  p 
is  the  local  thermodynamic  pressure;  1  is  the  unit  tensor;  /3 
is  the  second  coefficient  (or  bulk)  viscosity;  p  is  the 
first  coefficient  (or  shear)  viscosity;  and  5  is  the  strain 
tensor:  S  =  (Vv  +  (Vv)t)/2,  where  t  indicates  the  transpose, 
and  V  =  I,  e j 9/9x j ,  ej  being  a  cont ravar iant  (or  dual)  base 
vector  in  the  i ' th  direction.  (Note  that  the  gradient  of  a 
vector  is  a  tensor.)  For  an  ideal  gas  /3=0  and  for  most 
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1 1 

fluids  j3=0  [1,  p.  62],  so  that  equation  (2.1),  with  B=0  and 
p=constant,  can  be  simplified  to  the  well  known  form: 

pdv/dt  =  -Vp  +  pV2v  +  pV(V*v)/3  (2.2) 

To  ensure  a  total  energy  balance  throughout  the  fluid, 
it  is  necessary  to  satisfy  the  First  Law  of  Thermodynamics. 
This  law  states  that  for  an  elementary  volume  (AV)  as  it 
flows  along  its  path: 
dEt/dt  =  Q'  -  W' 

where  dEt  is  the  increase  in  total  energy  of  the  volume  in  a 
time  dt ,  Q'  is  the  rate  of  heat  transfer  to  the  volume,  and 
W ’  is  the  rate  of  work  done  by  the  volume  [1,  p.  265]. 

Here,  Q'  and  W'  are  not  written  as  dQ/dt  and  dW/dt,  because 
this  would  indicate  that  Q  and  W  could  be  transported  by  the 
fluid,  which  is  not  the  case.  Et  includes  kinetic  energy 
(pv*v/2)  and  internal  energy  (e).  Written  for  an  elementary 
volume:  Et  =  pAV(e  +  v*v/2).  Potential  energy  terms  are  not 
included  since  they  also  appear  in  W'  and  so  will  cancel  out 
[9,  p.  72-74].  Letting  Et=petAV,  with  et  being  the  total 
energy  per  unit  mass,  then 

pAVdet/dt  =  pAV(de/dt  +  d(v*v)/dt/2)  [1,  p.  266] 

Using  Fourier's  law  for  heat  flux,  the  rate  of  heat  flow 
into  the  elementary  volume  is  [9,  p.  73]: 

Q’  =  [ V •  ( /cVT )  ] AV  (2.3) 

where  k  is  the  constant  of  heat  conductivity,  and  T  is  the 
temperature  of  the  fluid.  The  rate  of  work  done  by  the 
fluid  is  [ 9 ,  p.  74 ] : 

W  =  [ -V* ( v • a ) ] AV 
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W'  =  [V»(pv)-V*(v*7r)]AV  (2.4) 

where  o  =  -pi  +  i  so  that  v  is  the  viscous  shear  stress 
tensor  7r  =  (-2/3pV*v)I  +  2pS.  The  First  Law  thus  becomes 
pdet/dt  =  [kV2T  -  V*(pv)  +  V«(v*ff)]  (2.5) 

In  the  systems  considered  in  this  project  the  fluids 
were  ideal  gases,  so  that  the  ideal  gas  law  p=pRT  was 
obeyed.  Other  fluids  could  be  modelled  by  a  particle 
approach  providing  their  equation  of  state  (i.e.  a  relation¬ 
ship  between  p,e,  and  T)  is  known.  Some  other  form  of  the 
energy  equation  (2.5),  involving  enthalpy  for  example,  could 
be  used  instead;  however,  some  equation  of  state  is  still 
needed . 

For  an  ideal  gas  de  =  Cv(dT),  where  Cv  is  the  specific 
heat  at  constant  volume.  Equation  (2.5)  can  thus  be 
written : 

dT/dt  =  [/cV2T  -  V-(pv)  +  V*  (v  •  7r )  ]/pCv 

-  (y*dv/dt)/Cv  (2.6) 

The  term  V»(y«j)  includes  the  effects  of  viscous  work 
v*(V*7r)  and  viscous  dissipation  0=7r :  Vv .  The  standard  form 
[1,  p.  267] 

dT/dt  =  (kV2T  -  pV*v  +  0)/pCv 

was  not  used  because  the  pressure  work  term  pV*v  cannot  be 
finite  differenced  conservatively,  and  led  to  large  scale 
instabilities  in  the  flow.  It  was  also  found  that,  even 
though  the  temperature  was  expected  to  be  constant,  equation 
(2.6)  was  necessary  for  the  channel  flow  simulation, 
otherwise  large  nonphysical  oscillations  developed  in  the 
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flow.  The  discharge  simulation  described  in  chapter  5  was 
stable  without  the  temperature  equation,  though  the  tempera¬ 
ture  was  not  expected  to  be  constant  there.  Note  in 
equation  (2.6)  that  the  term  v*dv/dt  is  particularly  suited 
to  the  particle  method  since  dv/dt  is  known  from  equation 
(2.2).  Note  also  that  a  mass  conservation  equation  is  not 
used  with  the  particle  approach  since  particles  are  not 
created  or  destroyed  in  a  closed  system,  so  that  mass  is 
naturally  conserved. 

2.2  Computational  model 

The  Lagrangian  equations  of  fluid  motion  for  an  ideal 
gas  (equations  [2.2]  and  [2.6])  can  be  modelled  by 
considering  the  fluid  as  composed  of  many  macroparticles, 
each  with  a  spatial  density  profile.  It  is  not  correct  to 
view  the  macroparticles  as  actual  physical  sections  of  the 
fluid,  since  the  particles  can  only  be  properly  described  by 
the  equations  which  relate  the  particle  quantities  to  the 
fluid  quantities  and  fluid  equations. 

In  the  systems  considered  herein,  each  particle  has 
associated  with  it  a  position,  a  velocity,  a  mass  density, 
and  a  temperature.  To  obtain  the  fluid  velocity,  tempera¬ 
ture,  and  density,  an  Eulerian  grid  is  used  on  which  the 
fluid  quantities  are  to  be  defined.  The  fluid  temperature 
and  velocity  at  the  middle  of  a  cell  (i.e.  at  a  grid  point) 
are  very  simply  defined  as  just  an  arithmetic  average  of  the 
temperatures  and  velocities  of  the  particles  inside  that 
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cell.  To  obtain  the  fluid  density  is  more  difficult  because 
of  the  assumed  finite  extent  of  the  particle  densities.  In 
this  work  the  particles  were  given  a  Gauss ian-shaped  density 
profile  (see  fig.  1.1).  It  would  be  possible  to  sum  the 
individual  density  contributions  of  all  the  macropart icles 
to  obtain  the  fluid  density  at  a  point,  but  this  would  be 
very  time  consuming.  Instead,  a  nearest  grid  point  method 
(NGP)  approximation  is  applied.  Defining  n ,  j  as  the  number 
of  particles  in  cell  i,j  whose  center  is  at  r,  ,,  then  the 
NGP  approximation  considers  all  particles  in  cell  i,j  as 
centered  at  £i ,  j .  The  density  at  any  other  other  grid  point 
m,n,  located  at  r_ m  n  ,  is  then  defined  to  be: 

p(rm  n,t)  =  p0 [ 1 1  ,  j  f (£m  i  .  j )n s  .  j (r i  jrt)]/ 

[I ,  ,  j  f (rm  n-r  j  j  )n0  ]  (2.7) 

where  pc  is  the  initial  fluid  density  associated  with  n0 
particles  per  cell;  the  sums  I,  }  are  over  all  grid  points; 
and  f ( r  -  £i,j)  describes  the  spatial  extent  of  a 
macropart icle  density  profile.  For  <Rn  ,  or  for  an 
n-dimensional  periodic  system,  the  denominator  of  equation 
(2.7)  is  written  as  2rn  2ann0  [6].  For  a  finite  system, 
however,  it  is  necessary  to  use  the  normalization  factor 
given  in  the  denominator  of  equation  (2.7)  since  this 
accounts  for  the  fact  that  near  boundaries  the  numerator  is 
less  in  value.  To  obtain  a  smooth  fluid  picture  from  the 
particles,  a  Gaussian  profile  of  width  ’a'  was  used  for 
f  (jr  -  £  j  j  )  i  .  e  . 

f ( rm  ,  n-r i  j )  =  exp[- (xm , n-x i  ,  j ) 2/( 2a 2 ) 

-  (ym.  n-y  i  ,  j )  V(2a2 )  ] 


(2.8) 
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where  'a2'  was  set  equal  to  between  one-half  and  one  times 
A2,  where  A  is  the  grid  spacing. 

The  fluid  quantities  are  thus  known  at  each  grid  point 
so  that  finite  differencing  can  be  used  to  evaluate  the 
right  hand  side  of  equations  (2.2)  and  (2.6),  and  approxi¬ 
mate  the  Lagrangian  derivatives  dv/dt  and  dT/dt  at  each  grid 
point  i,j.  As  each  particle  k  is  pushed  a  distance 
Axk  =  vkAt,  its  velocity  and  temperature  are  then  changed  by 
Avk  =  ( dv/dt ) i ,  j At ,  ATk  =  (dT/dt),,jAt  where  i,j  indicates 
the  NGP  to  the  particle.  With  this  new  configuration  of 
particles  the  fluid  quantities  can  be  recalculated  and  the 
particles  pushed  again,  etcetera. 

It  is  important  to  note  that  the  particles  do  not 
physically  extend  in  space,  since  there  are  two  major 
difficulties  in  such  an  interpretation.  First,  different 
parts  of  an  actual  spatial  particle  would  be  pushed  or 
heated  by  different  amounts  and  rates,  corresponding  to  the 
particle's  mass  in  different  cells  so  that  distortion  would 
occur  very  quickly  just  as  in  the  Lagrangian  finite  differ¬ 
ence  codes.  Secondly,  there  would  be  an  inconsistency  in 
assigning  the  fluid  quantities  (dv/dt)jfj  and  (dT/dt),  ti  to 
the  particles  in  cell  i,j.  This  is  because  these 
derivatives  are  based  on  the  spatial  derivatives  of  the 
fluid,  so  that  assigning  these  to  a  particle  implies  that  a 
truly  finite  particle's  velocity  and  temperature  would  vary 
in  space.  This  would  add  a  considerable  amount  of 
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complexity  to  the  problem.  The  particles  are  really  then 
just  a  means  of  carrying  the  Lagrangian  derivatives  on  an 
Eulerian  grid. 


2.2.1  The  particle  equations 

The  actual  form  of  the  equations  for  the  change  in  the 
kth  particle's  position,  velocity  and  temperature  are: 

dXk/dt  =  vk  (2.9) 

dvk/dt  =  [L(v,T,p) j  f j ]/p i , j  (2.10) 

dTk/dt  =  [E(v,T,p) j  ( j]/Pi , jCv  -  ( v k • dv k /dt ) /Cv  (2.11) 
Li ( j  and  Ej  , j  indicate  the  values  of  the  right  hand  side  of 
equations  (2.2)  and  (2.5)  respectively,  finite  differenced 
at  the  cell  i,j  that  the  particle  is  in. 

All  quantities  in  L  and  E  should  be  staggered  by  a  half 
time  step  from  the  values  used  in  the  approximations  to 
dv/dt ,  dT/dt  [2,  p.  405].  A  correct  time-centered  scheme 
for  equations  (2.9)— (2. 1 1 )  written  at  t=nAt  is: 
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where  the  fluid  quantity  vv  "  '  is 
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(2.15) 


-  similarly  for  T(n).  Note  that  equation  (2.15)  is  valid 
only  for  At  being  a  constant.  The  sums  Lk ( j  ( j >  are  over  all 
particles  in  the  cell  i , j ,  and  (ni(j)n  =  (£k(i(j)1)n  is  the 
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number  of  particles  in  the  cell  i,j  at  time  t=nAt.  This 
peculiar  form  for  v(n)  can  be  written  as 

(v(n))j,j  =  (Zk ( s  , j ) v£ n } ) /( n j  ; j ) n  which  is  simply  a  linear 
approximation  to  the  average  of  the  particle  velocities  in 
cell  i,j  at  t=nAt.  The  form  of  (vk‘dvk/dt)n  in  equation 
(2.14)  uses  the  value  of  (dvk/dt)n  from  equation  (2.13).  An 
interpolation  is  used  here  for  the  kth  particle  velocity  vk  , 
rather  than  the  extrapolation  which  had  to  be  used  for  the 
fluid  velocity  (y(n))i,j  in  L  and  E.  This  is  because,  at 
the  point  in  the  code  where  equation  (2.14)  is  used  to  heat 
or  cool  the  particles,  vk+l/2  is  known  from  equation  (2.13). 
In  calculating  L  and  E  however,  v" + 1 '  2  was  not  known  yet  and 
so  extrapolation  for  v(n)  (2.15)  was  required.  In  (2.13) 
and  (2.14),  pn  is  the  value  of  the  density  based  on  the 
positions  of  the  particles  x£  at  t=nAt.  This  time-centering 
scheme  improves  the  accuracy  of  the  values  of  d/dt  by 
assuming  the  quantities  vary  linearly  between  time  steps. 

With  the  form  of  equations  ( 2  .  1  2 ) - ( 2 . 1 4 )  it  is  possible 
for  a  particle  velocity  or  temperature  to  be  quite  different 
from  the  fluid  velocity  or  temperature  of  the  cell  the 
particle  is  in.  This  is  referred  to  as  multistreaming,  and 
can  be  removed  by  adding  an  artificial  drag  term  [2]  to  the 
particle  equations  (2.13)  and  (2.14): 

yr  1/2  -  yr  1/2  =  [L"  (  j  ]  At/p"  ,  j  -  am  (v  k  ~ 1  ;  2  -<v> n )  (2.16) 
T  k  +  1 7  2  -  TE“1/2  =  [E? ,  j ]At/p? ,  jCv  -  (vk  *dvk/dt ) n/Cv 

-  a  t  (TJJ-  1  '  2-<T>n  )  (2.17) 

where 
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<v>"  =  (Ik(i, j)vk"-’/2)/(ni , j ) "  (2.18) 

(similarly  for  <T>n)  and  am  and  at  are  constants  (referred 
to  collectively  as  "a" )  0^a<1.  The  value  of  am  and  at  used 
must  be  examined  carefully,  as  discussed  in  section 
4.5.  <v>"  and  <T>n  appear  similar  to  v(n)  and  T(n)  from 

(2.15),  but  the  peculiar  mix  of  time  steps  in  (2.18)  is 
necessary  so  that,  when  summed  over  all  particles  in  a  cell, 
the  drag  terms  vanish.  That  is,  when  the  last  term  on  the 
right  hand  side  of  equation  (2.16)  is  summed  over  all 
particles  in  cell  i , j ,  the  result  is  -am(<v>n  -  <v>")  =  0. 
This  ensures  that  in  any  cell  no  net  force  appears  from  the 
artificial  drag  terms  [2]. 

2.2.2  Finite  differencing  of  the  equations 

Consider  now  the  way  in  which  to  approximate  the  terms 
appearing  on  the  right  hand  side  of  equations  (2.2)  and 
(2.5),  which  comprise  Lj , j  and  Ei(j  in  the  particle 
equations  (2.16)  and  (2.17).  The  fluid  is  divided  into 
cells  by  the  Eulerian  grid;  at  the  middle  of  each  cell  i,j, 
the  fluid  quantities  vi;Jn)  and  Ti(5n)  are  known,  and  these 
can  be  considered  as  the  average  value  of  the  fluid  quan¬ 
tities  for  that  cell.  The  method  chosen  to  approximate  the 
derivatives  should  be  conservative  between  time  steps  in  the 
total  momentum  and  energy  of  the  system.  That  is, 

Ikmkdvk/dt  (where  the  sum  is  over  all  particles)  must  be 
nonzero  only  due  to  momentum  being  added  or  lost  at  the 
boundaries,  assuming  there  are  no  body  forces.  If  there  are 


19 


body  forces,  then  this  statement  must  still  be  true  if  the 
contributions  of  the  body  force  B  to  the  sum  are  left  out. 
Similarly  in  energy,  Ikmk (de ,/dt ) k  must  be  nonzero  only  due 
to  energy  transfer  at  the  boundaries,  if  the  contributions 
from  any  internal  body  heat  sources  (for  example  Joule 
heating  in  an  MHD  fluid)  are  left  out  of  the  sum.  (Note: 
strictly  speaking  it  is  I j  ; j [ p3v/3 t ] j , j ,  and 
I i ( j [ pde t /3t ] j , j ;  where  the  sums  are  over  all  grid  points 
and  the  derivatives  are  partial,  which  should  be  conserved; 
but  this  is  not  possible  with  a  macroparticle  scheme.  By 
using  total  derivatives  in  the  sums,  the  system  will  be 
conservative  in  a  Lagrangian  rather  than  an  Eulerian  sense.) 
Using  equation  (2.16)  then: 

Ikmkdvk/dt  =  I j  , jZk ( j , j )L(v,T, p) i ( jmk/p j ( j  (2.19) 

The  sum  Zk  (  ■,  ;  }  }  is  over  all  particles  k  in  the  cell  i,j. 

The  dummy  index  k  appearing  on  the  right  hand  side  of  (2.19) 
has  no  connection  with  the  dummy  index  appearing  on  the  left 
hand  side.  The  mass  of  each  particle  mk  on  a  uniform  grid 
is  mk  =  pi(jA3/ni(j,  where  A  is  the  grid  spacing.  So  (2.19) 
becomes : 

Ikmkdvk/dt  =  A3Z j  , jL(v,T,p) j  , j  (2.20) 

The  energy  equation  can  be  treated  in  the  same  way  to  yield: 

Lkmk (de t/dt ) k  =A3Ii(jE(y,T,p)i(j  (2.21) 

From  equations  ( 2 . 20 )  ,  ( 2 . 2 1 )  it  is  apparent  that  we  need  a 
method  of  finite  differencing  which,  when  applied  to  the 
right  hand  side  of  equations  (2.2)  and  (2.5),  yields  sums 
IjjL  and  Ij(jE  that  are  dependent  only  on  the  values  of  the 
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fluid  quantities  at  the  boundaries.  In  cartesian 
coordinates  such  a  method  exists  and  will  be  termed  "flux 
conservative"  differencing  [3],  There  is  no  physical 
interpretation  to  the  name.  This  method  can  be  applied  to  L 
and  E,  but  to  apply  it  to  the  viscous  terms  in  E  is  tedious. 
The  term  V«(v*tt),  involving  viscous  work  and  viscous 
dissipation  in  equation  (2.5),  was  expected  to  be  small  in 
comparison  to  the  other  terms  and  was  not  of  prime  import¬ 
ance  in  the  geometries  considered;  therefore,  a 
nonconservative  form  was  used  for  it.  All  other  terms  in  E 
and  L  were  differenced  by  flux  conservative  methods.  In 
noncartesian  coordinates  it  is  not  possible  to  difference 
all  terms  flux  conservatively;  however,  a  "force 
conservative"  differencing  can  be  introduced  in  an  attempt 
to  have  a  stable  form  of  differencing. 

The  principle  of  flux  conservative  differencing  is  to 
have  the  values  of  L  and  E  cancelled  between  cells  in  the 
Eulerian  grid.  For  example,  the  term  9(FG)/3x1,  where  F  and 
G  represent  fluid  quantities  or  position,  can  be  differenced 
on  a  uniform  grid  as: 

[ 3 ( FG ) /3x i ] i  <  j  =  (Fj  +  i/  2,jGj  +  i/  2,j  -  Fj_i/2,jGj_i/2,j)/A 

where  Fi  +  1/2,j  =  ( F  ■,  *  i  ,  j  +F  ;  ,  i  )  /2  ,  and  i,j  refers  to  the 
position  of  the  grid  point  i,j  on  the  mesh.  This  is  easily 
seen  to  satisfy  the  form  of  flux  conservative  differencing 
when  summed  over  the  entire  mesh.  When  G  is  a  derivative: 

G  =  3H/3x1f  this  becomes: 

[3(F3H/3x, )/8x, ] , , ,  =  [F , . , , , , j (H , . , , j -H , , j ) 

-  F, . ,/2 , ,(H, , j-H, . , . j) ]/A* 
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Some  other  flux  conservative  finite  difference 
approximations  are  [3,  p.  59]: 

OF/3X,  )  ,  ,  j  =  (F,  «,  ,  ,  j-F,  .  ,  ,  j  )/2A 

(3:F/3xj)i,j  =  (Fi  +  ,  ,  j -2F  ;  ,  j+F,  . 1  ,  j )/A2 

[3(9G/ax1  )/Sx2]  ,  ,  jS[  OG/0X,  )  |  i  ,  j  *  ,-(3G/3x,)  |  ,  ,  j.  ,  ]/2A 

—  ( G i ♦ i ,  j ♦ i “G  j  _ i  ,  j  +  i  ~G  j  +  i  ,  j-i +G i_i(j_i  )/4A2 

2.2.2. 1  Modified  qradp  differencing 
In  differencing  the  Vp  term  in  equation  (2.2),  ordinary 
flux  conservative  differencing  proved  unsatisfactory  because 
it  led  to  large  scale  clumping  of  the  particles,  producing 
large  variations  in  density  over  only  a  few  grids.  Particle 
oscillations  were  created  by  these  density  fluctuations  and 
large  scale  random  velocities  resulted;  this  is  the  source 
of  the  very  undesirable  noise  levels  found  in  previous  NGP 
codes.  As  a  simple  example  of  this,  consider  a  cell  which 
has  a  large  number  of  particles  and  thus  a  high  mass 
density.  Suppose  this  cell  is  surrounded  by  cells  which  all 
have  the  same  number  of  particles,  and  this  number  is  less 
than  that  of  the  middle  cell.  Then  the  particles  in  the 
middle  cell  will  erroneously  see  no  Vp  force  since: 

<Pi  , j ♦ i“Pi  , j - i )/2A  =  (Pi ♦ i , j  — p i - i , j )/2A  =  0;  this  is  clearly 
in  error  since  in  fact  the  particles  in  this  cell  should 
feel  a  Vp  force  which  pushes  them  away  from  the  middle  of 
the  cell  in  order  to  fill  in  the  lower  density  region  sur¬ 
rounding  them. 
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The  difficulty  introduced  by  a  conventional  finite 
differencing  of  Vp  evaluated  directly  at  the  cell  center  was 
removed  by  evaluating  Vp  at  the  midpoints  between 
neighboring  cell  centers,  and  then  linearly  interpolating 
between  these  values  based  on  a  particle's  position  in  the 
cell.  For  example,  in  cartesian  coordinates  x  and  y,  for 
the  kth  particle  in  cell  i,j  (see  figure  2.1): 

Op/dx)  k  =  Op/dx  |  i  -  i  /  2  ,  j 

+  ( 9p/9x |  j  +  , / 2 , j-9p/9x |  j _  !  /  2  ,  j ) (xk-x j _ , / 2 ,  j )/A  (2.22) 
where  9p/3 x |  i _ i / 2 ,  j  =  (pi , j”Pi - i , j )/A,  and  similarly  for  the 
other  derivative  term.  (xk-Xj-i/2, j)  is  the  distance  to  the 
right  of  xi_1/2,j  ,  the  midpoint  between  cell  centers  i  —  1 , j 
and  i , j .  Op/Oy  is  evaluated  in  a  similar  manner  between 
midpoints  i,j+1/2  and  i ,j— 1/2.  All  forces  could  be  assigned 
to  the  particles  in  this  way,  however  it  was  found  suf¬ 
ficient  to  deal  with  just  the  Vp  term  in  this  way. 

Note  that  equation  (2.20)  can  no  longer  be  strictly 
satisfied,  since  the  value  of  Vp  in  L  for  each  particle  in 
cell  i , j  is  no  longer  the  same.  Equation  (2.20)  must  now  be 
written : 

Ikmkdvk/dt  =  A3I i , jlk ( i , j ) [ (-Vp) k/n j  , j ]  +  A3Ii(j(L  +  Vp) 
The  second  sum  on  the  right  hand  side  vanishes  since  terms 
other  than  Vp  in  L  are  differenced  flux  conservatively.  The 
first  sum  however  is  not  identically  zero  unless  Vp  is  the 
same  for  all  ni(j  particles  in  the  cell  i,j.  This  should 
not  cause  any  instabilities  to  form  however  since  when 
averaged  over  many  time  steps  the  sum  will  be  approximately 
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Figure  2.1  Vp  finite  differencing.  9p/9x  is  evaluated  at  the 
crosses(x)  and  9p/9x  inside  cell  i,j  is  then  a 
linear  interpolation  in  x  between  these  two 
values  as  in  equation  (2.22).  9p/9y  is  evaluated 
at  the  circles(O)  and  is  similarly  interpolated. 
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zero  for  steady  state  conditions.  This  is  because  the 
individual  (Vp)k  are  obtained  by  flux  conservative  means. 

To  show  that  the  sum  is  approximately  zero,  consider  for 
simplicity  just  the  x  direction  in  cartesian  coordinates. 
Then  the  particle  k  from  equation  (2.22)  contributes 
Op/9x)k  when  it  occupies  position  xk.  Assuming  steady 
state  conditions,  then  if  at  some  other  time  a  particle 
occupies  the  position  xm=-xk,  and  at  some  other  time  another 
particle  occupies  the  position  xs=-xk+2A,  then  the  contri¬ 
bution  to  the  Vp  sum  from  particle  k  will  be  zero.  Because 
there  are  many  particles  per  cell  these  events  are  likely  to 
be  approximately  true  after  relatively  few  time  steps. 

2. 2. 2. 2  Force  conservative  differencing 
As  mentioned  above,  some  terms  in  L  and  E  cannot  be 
differenced  flux  conservatively  when  L  and  E  are  written  in 
noncartesian  coordinates.  All  such  terms  are  of  the  form 
F9G/9x,  and  cannot  be  differenced  so  that  their  sum  over  the 
grid  vanishes.  The  value  of  F9G/9x  is  instead  approximated 
at  a  cell  center  by  an  average  of  a  finite  difference 
approximation  at  the  midpoints  between  cell  centers,  i.e.  a 
difference  of  the  form 

F9G/9x  =  [Fj+1/2,j(Gj+1)j ~G i ,  j ) /A  +  Fj_i/2lj(Gj(j_Gj-i(j 
is  used.  This  is  done  in  order  to  have  the  values  of  F9G/3x 
continuous  at  the  interfaces  of  the  cells.  The  term  "force” 
has  no  physical  significance  in  the  expression  "force 
differencing"  . 


)/A]/2 
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Note  that  a  term  of  the  form  3 (F9G/9x , ) /9x 2  is  flux 
conservative,  even  though  the  term  F9G/9X!  is  differenced  as 
a  force  conservative  term,  i.e. 

9(F9G/9x1)/9x2s{[Fi+1/2(Gi+1-Gi )/A  +  Fi_1/2(Gi-G,_1)/A]j+1/2 
~  t  F j  +  i / 2 (G i  +  i ~G i ) /A  +  F  j  _ i / 2 (G  j -G i-i)/A]j_1/2}/ 2 A 

where 


i  +  1  /  2  ,  j  +  1  /  2 


=  [(F, 


+  1  ,  j  +  1 


+F 


)/2  +  (F 


i  +  1 


+F i , j )/2]/2 


2.2.3  Variable  grids 

One  can  apply  the  above  method  of  finite  differencing 
to  a  variable  grid  [3,  p.  60];  however,  the  macroparticle 
method  becomes  more  expensive  in  this  case,  not  only  because 
of  an  increase  in  the  complexity  of  the  finite  differencing 
and  calculating  the  density  on  a  variable  grid,  but  also 
because  there  is  a  minimum  number  of  particles  per  cell 
required  to  maintain  a  reasonable  accuracy.  Since  the 
number  of  particles  per  cell  is  proportional  to  a  cell  area, 
then  the  larger  cells  require  more  particles  per  cell  than 
the  minimum.  This  adds  to  the  computational  work  necessary 
since  the  larger  cells  have  more  particles  than  is  necessary 
for  good  accuracy.  Using  a  variable  grid  to  obtain  better 
resolution  in  certain  areas  may  thus  not  be  as  cost 
effective  as  in  conventional  fluid  codes.  This  is 
especially  true  if  most  cells  are  of  a  size  significantly 
larger  than  the  smallest  cell.  A  uniform  grid  was  used  for 
all  problems  considered  in  this  project  since  there  was 
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little  need  for  a  variable  grid  for  the  geometries  under 
consideration . 

2.2.4  General  conservative  form  of  the  equations 

In  order  to  be  able  to  satisfy  the  condition  that  the 
change  in  total  momentum  and  total  energy  depend  only  on 
fluid  quantities  at  the  boundaries,  the  form  of  the  momentum 
equation  to  be  finite  differenced  (2.2)  must  be  of  the 
general  conservative  form  (ignoring  body  forces  as  mentioned 
earlier ) : 

pdv/dt  =  div(fi)  (2.23) 

where  div (M)  =  M* V  and  the  derivatives  operate  on  M .  Note 
that  for  a  symmetric  tensor  div(fi)  =  V*M  and  all  tensors  in 
the  previous  equations  are  symmetric.  Integrated  over  the 
entire  volume  (V)  this  equation  yields: 

/( pdv/dt )dV  =  /div(fi)dV  =  Jft*hdS  (2.24) 

where  S  is  the  surface  enclosing  the  system  and  n  is  the 
unit  normal  to  the  surface.  In  this  way,  the  only  physical 
sources  or  sinks  of  pdv/dt  result  from  the  system 
boundaries.  Similarly,  the  energy  equation  (2.5)  must  be  of 
the  form  (ignoring  internal  heat  sources  as  mentioned 
earlier ) : 


pdet/dt  =  V • A 

(2.25) 

where  A  is  a  vector,  so  that 

/ ( pde  t /dt ) dV  =  / A • hdS 

(2.26) 

Although  in  noncartesian  coordinates  the  sums  on  the  right 
hand  side  of  (2.20)  and  (2.21)  cannot  be  made  dependent  only 
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on  the  boundary  conditions,  it  is  still  desirable  to  have 
the  equations  in  conservative  form  ( [ 2 . 23 ] , [ 2 . 25 ] )  so  that 
the  integral  equivalent  of  ( 2 . 20 )  ,  ( 2 . 2 1 )  ( [ 2 . 24 ] , [ 2 . 26  ] ) 
results  in  there  being  no  nonphysical  sources  or  sinks  of 
momentum  or  energy. 

2.2.5  Boundary  conditions 

The  boundary  conditions  applied  to  the  finite 
differenced  form  of  L  and  E  are  the  "ghost  point"  or  "guard 
point"  conditions.  That  is,  the  boundaries  of  the  fluid 
system  are  made  to  lie  midway  between  two  grid  points,  one 
interior  and  the  other  exterior  to  the  system,  the  exterior 
grid  points  being  where  the  boundary  conditions  are  defined. 
The  values  of  the  fluid  quantities  at  the  fictitious 
exterior  grid  points  are  defined  so  that  the  correct 
boundary  conditions  for  the  problem  being  considered  are 
satisfied.  For  example,  if  the  fluid  velocity  is  to  be  zero 
on  a  containing  wall,  then  the  ghost  point  fluid  quantities 
would  be  the  negative  of  those  interior  to  the  wall.  To 
have  a  Dirichlet  condition  at  a  boundary,  then  the  value  of 
the  quantity  at  an  exterior  grid  point  is  set  so  that  the 
average  of  the  ghost  point  value  and  the  value  at  the  grid 
point  interior  to  the  boundary  is  the  Dirichlet  value.  To 
have  a  zero  derivative  of  some  quantity  at  a  boundary,  then 
the  ghost  point  value  is  made  equal  to  the  interior  value. 
All  boundary  conditions  for  the  fluid  equations  can  be 
satisfied  in  this  manner.  In  addition,  at  a  rigid  wall  the 
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particles  are  elastically  reflected. 
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CHAPTER  III 


Experimental  and  Analytical  Channel  Flow 

In  developing  the  viscous  fluid  particle  method,  a  very  well 
studied  geometry  was  needed  for  simulation  so  that  compari¬ 
son  of  code  results  to  experimental  results  could  be  done 
reliably.  For  this  reason,  the  flow  of  a  fluid  between  two 
very  large  smooth  planar  plates,  known  as  smooth  channel 
flow,  was  selected.  The  quantity  of  available  experimental 
data  and  analytical  results  for  channel  flow  is  large,  and 
only  a  small  amount  will  be  discussed  here.  In  verifying 
the  code  results,  several  important  parameters  which 
characterize  the  flow  were  chosen  for  comparision,  although 
many  more  less  critical  parameters  could  also  be  compared. 

In  describing  the  flow  of  any  fluid  it  is  important  to 
distinguish  between  laminar  flow,  which  is  free  from  any 
turbulence,  and  turbulent  flow.  According  to  Hinze  [11,  p. 
2],  "Turbulent  fluid  motion  is  an  irregular  condition  of 
flow  in  which  the  various  quantities  show  a  random  variation 
with  time  and  space  coordinates,  so  that  statistically 
distinct  average  values  can  be  discerned."  Most  realistic 
flows  are  turbulent  and  in  the  following  discussion  only 
turbulent  flows  will  be  considered. 

3.1  Experimental  and  analytical  relations  for  a  few 

variables 

Since  there  is  a  random  component  to  turbulence,  a 
turbulent  flow  can  only  be  described  consistently  in  terms 
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of  time  averaged  quantities.  It  was  necesary  to  time 
average  the  code  results  to  compare  to  other  results  since 
the  instantaneous  fluid  equations  were  simulated.  Note  that 
by  time  averaging  the  fluid  equations  (2.2)  and  (2.6),  the 
so  called  Reynolds  averaged  equations  for  the  time  average 
fluid  quantities  can  be  obtained  [11,  p.  21].  These 
equations  contain  too  many  variables  to  be  solved  without 
some  assumptions  being  made,  and  are  not  suitable  for  a 
particle  simulation  because  the  particles  behave  in  an 
instantaneous,  not  a  time  averaged  manner.  The  time 
averaged  quantities  can  always  be  obtained  from  the 
instantaneous  values  by  a  straight-forward  time  averaging. 

In  most  literature  a  time  averaged  quantity  is  indicated  by 
a  line  above  the  instantaneous  quantity.  Unfortunately,  due 
to  word  processing  limitations  this  notation  could  not  be 
used.  Instead,  in  the  following  discussion,  the  triangular 
brackets  <  >  enclosing  a  quantity  will  indicate  the  time 
averaged  value  of  that  quantity  unless  otherwise  specified. 
The  time  average  is  assumed  to  be  taken  for  a  long  enough 
period  that  the  time  averaged  quantity  is  no  longer  a 
function  of  time.  A  capital  letter  U  will  be  used  to 
indicate  <vx>,  the  time  averaged  along  stream  velocity 
component;  and  a  capital  letter  P  will  be  used  to  indicate 
<p>,  the  time  averaged  thermodynamic  pressure.  Other 
conventions  will  be  described  as  they  are  encountered. 

For  channel  flow  a  typical  time  averaged  along  stream 
velocity  profile  is  shown  in  figure  3.1.  Note  that  the 
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Figure  3.1  A  typical  profile  of  along  stream  velocity  U  as  a 
function  of  cross  stream  distance  y  for 
turbulent  flow  through  a  channel  of  cross  stream 
width  D. 
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profile  is  symmetric  about  its  centerline.  Channels  of  dif¬ 
ferent  width,  different  midstream  velocity,  or  different 
fluids  will  have  different  profiles;  however,  it  is  known 
that  if  the  quantity  (U0_U)/u*  is  plotted  against  cross 
stream  distance  2y/D  where  D  is  the  width  of  the  channel, 
then  a  profile  which  is  nearly  universal  for  all  turbulent 
channel  flows  is  obtained  [1,  p.  587].  u*  is  called  the 
wall-friction  velocity  and  is  a  characteristic  constant  of 
the  flow  defined  by: 

(u* ) 2=t 0/p  (3.1) 

where  r0  is  the  value  of  the  time  averaged  x-direction 
shearing  stress  felt  at  the  wall;  U0  is  the  time  averaged 
midstream  velocity;  and  U  is  the  time  averaged  along  stream 
velocity.  This  profile  is  somewhat  different  than  that 
obtained  for  a  flat  plate,  but  is  very  similar  to  that  for  a 
round  pipe.  When  very  near  the  wall,  the  velocity  profiles 
for  the  flat  plate,  round  pipe,  and  channel  are  the  same 
[11,  p.  717].  Analytic  functions  for  (U0“U)/u*  can  be 
obtained  [1,  p.  586], [12],  but  these  are  always  verified  by 
comparison  to  experimental  results  so  that  experimental 
results  were  used  for  comparison  of  the  code  results.  The 
data  of  Nikuradse  and  Donch  obtained  from  [12,  p.  482]  was 
used  and  is  plotted  in  figure  3.2. 

It  is  also  possible  to  define  velocity  profiles  U/u* 
which  are  functions  of  yu*/v ,  where  u=n/p,  but  these  are 
only  valid  for  a  fraction  of  the  channel  width.  Another 
analytic  function  which  is  somewhat  approximate,  but  is 
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Figure  3.2  Experimental  plot  of  universal  velocity  profile 

(U0~U)/u*  versus  normalized  distance  from  the 
wall?  obtained  from  [12,  p.  482] 
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often  used  because  of  its  simplicity,  is  the  power  law  [7, 
p.  723] 

U/U0  =  (2y/D) 1 /n  (3.2) 

where  n  is  a  parameter  which  varies  with  the  Reynolds  number 
(based  on  the  midstream  velocity)  Re0=pDU0/M.  The  Reynolds 
number  represents  the  ratio  of  inertial  forces  to  viscous 
forces  [1,  p.  13],  and  flows  which  have  the  same  Reynolds 
number  are  dynamically  similar. 

In  order  to  compare  the  code  results  to  the  data  of 
Donch  and  Nikuradse  (fig.  3.2),  it  is  necessary  to  obtain 
the  value  of  u*  for  the  flow  conditions  being  simulated. 

This  can  be  done  by  using  the  definition  of  the  resistance 
coefficient  of  pipe  flow  [1,  p.  612]: 

AP/L  =  Xp ( <U>  s ) 2 /2h  (3.3) 

where  X  is  the  resistance  coefficient;  AP  is  the  drop  in  P, 
the  time  averaged  pressure,  over  an  along  stream  distance  L; 
<U> s  is  the  value  of  U  averaged  across  the  channel;  and  h  is 
the  hydraulic  diameter  h=4A/C,  A  denoting  the 
cross-sectional  area  and  C  denoting  the  wetted  perimeter 
-  for  channel  flow  h=2D. 

Now,  the  forces  acting  on  a  section  of  fluid  of  length 
L,  width  D,  and  of  unit  depth  must  be  zero  for  equilibrium 
flow,  so  that  the  total  shear  force  due  to  r0  acting  on  the 
fluid  at  the  wall  must  equal  the  total  pressure  force  acting 
over  the  cross-section  due  to  the  pressure  drop  AP.  Note 
that  P,  the  time  averaged  pressure,  cannot  vary  over  the 
cross  section,  for  obvious  physical  reasons.  In  the  form  of 
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an  equation  this  statement  of  zero  total  force  can  be 
written : 

DAP  =  2Lr0  (3.4) 

By  combining  (3.1),  (3.3)  and  (3.4)  then: 

u*  =  <U>SX’ 7  2//8  (3.5) 

X  is  a  function  of  the  Reynolds  number  Red=Dp<U>s/p  and  can 
be  obtained  graphically  from  a  "Moody  diagram",  which  is 
essentially  an  experimental  plot  of  X  versus  Red;  or  an 
analytic  approximation  to  the  Moody  diagram  can  be  used. 

The  latter  route  was  chosen  and  X  was  approximated  as  [13]: 

X  =  [  0 . 868-^n  (Red/(  1 . 964^nRed  -  3.82))]‘2  (3.6) 

for  smooth  channel  flow.  In  this  way,  since  <U>S  is  easily 
obtained  from  the  code  results  of  U(y),  the  value  of  u*  is 
known  so  that  the  code  results  of  U  can  be  compared  to  the 
experimental  results  (U0-U)/u*. 

From  equation  (3.3),  it  is  also  clear  that  the  pressure 
P  is  an  exact  linear  function  of  downstream  distance  L  for 
an  incompressible  flow  ( p^constant ) .  For  a  compressible 
flow  for  which  <U>S  and  p  do  not  vary  appreciably  over  L 
then  equation  (3.3)  is  still  approximately  true,  and  (3.3) 
can  be  written: 

P  =  LXp<U> 2 /4D  +  constant  (3.6A) 

By  applying  the  same  argument  which  yielded  equation 
(3.4),  but  now  to  a  section  of  fluid  centered  on  the  channel 
centerline  which  is  of  arbitrary  width  (rather  than  of  width 
D),  then  it  is  clear  that  the  time  averaged  shear  stress  r 
is  a  linear  function  of  distance  from  the  centerline.  From 


36 


this  it  is  also  found  that  r  is  also  symmetrical  about  the 
centerline  and  is  zero  on  the  centerline.  Written  for  the 
lower  half  channel  then  [7,  p.  717]: 

r/r o  =  1  -  2y/D  (3.7) 

where  y  is  the  distance  from  the  lower  wall. 

If  the  fluid  equations  are  written  in  cartesian 
coordinates  and  time  averaged  as  mentioned  earlier,  then  the 
shear  stress  r  appears  in  the  boundary  layer  approximation 
to  the  Reynolds  averaged  x-momentum  equation  as  [1,  p.  563]: 

p(U9U/9x+V9U/9y )  =  -9P/9x  +  9r/9y  (3.8) 

where  r  has  the  form  [1,  p.  563]: 

r  =  p9U/9x  -  p<uv>  (3.9) 

u  and  v  are  the  instantaneous  turbulent  fluctuations  of  the 
total  instantaneous  velocity  about  the  time  averaged 
velocity,  i.e.  vx=U+u,  vy=V+v,  (vz=W+w).  The  quantity  <uv> 
is  the  time  averaged  value  of  uv,  and  (-p<uv>)  is  the 
familiar  turbulent  Reynolds  shear  stress  [1,  p.  563]. 
(-p<uv>)  is  not  a  true  stress  but  is  rather  a  momentum 
exchange  effect.  That  <uv>  is  not  zero  is  a  well  known 
phenomena  of  turbulent  shear  layers.  A  somewhat 
oversimplified  explanation  of  this  fact  can  be  given  by 
considering  figure  3.3.  It  is  seen  that  if  a  portion  of  the 
fluid  in  a  shear  layer  moves  to  a  region  where  U  is  higher, 
as  a  result  v>0,  then  once  it  arrives  at  its  new  position  it 
will  have  created  a  u<0  in  the  flow  there.  Here  the 
oversimply ing  assumption  must  be  made  that  the  velocity  of 
the  section  of  fluid  that  moved  upward  was  preserved  as  it 
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Figure  3.3  In  a  turbulent  shear  flow  along  a  wall  <uv>  is 

less  than  zero.  See  page  36  for  an  explanation. 
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moved  the  small  cross  stream  distance  to  its  new  position. 
For  v<0  then  u>0  will  be  created,  so  that  on  average  uv  is 
less  than  zero.  The  viscous  part  of  r  in  equation  (3.9)  is 
the  usual  viscous  shear  stress  associated  with  molecular 
viscosity  effects,  and  is  experimentally  found  to  be 
negligible  in  comparison  with  the  Reynolds  shear  stress 
(-p< uv>)  when  beyond  a  short  distance  from  the  wall. 

Since  by  definition  <u>=<v>=0,  then  it  is  necessary  to 
define  the  quantities  <v2>1/2/U0  and  <u2>1/2/U0  to  obtain  a 
measure  of  the  intensity  of  turbulence.  The  intensity  of 
turbulence  varies  across  the  channel,  however  a  typical 
value  of  <u2>172/U0  for  channel  flow  is  5-10%.  <v2>1/2/U0 

is  somewhat  less  than  this,  being  about  75%  of  the 
x-direction  value  [11,‘p.  639,  p.  725].  A  better  estimate 
for  <u2>  can  be  obtained  using  figure  7-60  from  [11],  from 
which  one  can  estimate  that  <q2>/(u*)2  =  4,  where 
<q2>=<u2>+<v2>+<w2>;  but,  <v2>=<w2>=0 . 75<u2> ,  so  this  yields 
<u2>  =  1.6(u*)2.  Substituting  in  equation  (3.5)  for  u* , 
then  we  have: 

<u2> 1 7  2/<U>  s  s  1 . 26 ( X/8 ) 1 7  2  (3.9A) 

3.2  Definition  of  layers  in  channel  flow 

For  convenience,  the  channel  is  usually  divided  into 
several  regions  based  on  distance  from  the  wall.  The  region 
closest  to  the  wall,  where  the  flow  is  directly  affected  by 
the  condition  at  the  wall,  is  called  the  ’’inner"  or  "wall" 


39 


region  [11,  p.  588].  Beyond  this  region  towards  the  center, 
the  flow  is  only  indirectly  affected  by  the  wall  and  is 
called  the  core  region.  Very  close  to  a  perfectly  smooth 
wall  there  is  also  a  "viscous  sublayer"  and  a  "buffer" 
layer.  The  viscous  sublayer  is  a  very  thin  layer  where  the 
flow  is  predominantly  viscous,  between  the  wall  itself  and 
the  inner  region.  The  buffer  region  is  the  name  given  to 
the  transition  region  between  the  viscous  sublayer  and  the 
inner  region.  It  can  be  shown  [11,  p.  718,  p.  615]  that  the 
velocity  profile  is  a  linear  function  of  y  in  the  viscous 
sublayer,  U  being  zero  for  y=0.  It  is  also  known  that  in 
the  viscous  sublayer  p3U/3y>> ( -p<uv> ) ,  so  r  is  essentially 
due  to  viscosity  here.  In  contrast,  in  the  core  region  the 
Reynolds  stress  is  very  much  larger  than  the  viscous  shear 
stress.  Throughout  the  buffer  and  inner  layers,  r  is  due  to 
both  the  viscous  and  turbulent  shear  stresses. 

Defining  the  normalized  cross  stream  distance  y+=yu*/v, 
then  the  viscous  sublayer  generally  extends  from  the  wall  to 
y+  ~  5.  In  the  case  of  flow  past  a  flat  plate,  or  high 
Reynolds  number  pipe  and  channel  flow,  the  inner  and  buffer 
regions  intersect  at  y+  ~  30  and  the  inner  wall  and  outer 
regions  usually  overlap  near  y+  ~  1000;  however,  for  low 
Reynolds  number  flow  the  channel  half-width  may  be  less  than 
these  distances.  Since  the  inner  region  usually  only 
occupies  one-tenth  or  so  of  the  channel,  then  in  this  case 
the  distance  y*  from  the  wall  to  the  boundaries  between 
these  regions  must  depend  on  the  flow  parameters. 
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3.3  Eddy  size  analysis 

In  modelling  turbulent  fluid  flow  on  a  computational 
grid,  it  is  necessary  to  have  the  cell  sizes  about  the  same 
size  as  the  smallest  eddies  of  interest.  We  chose  to 
simulate  most  of  the  range  of  eddy  sizes  in  the  channel,  so 
that  the  approximate  size  of  the  smaller  eddies  was 
required.  This  size  is  known  to  be  the  Taylor  microscale 
Xz  ,  defined  by  (Xz ) 2=2<u2>/<  Ou/9y ) 2>.  Note  that  it  is  con¬ 
ventional  in  the  literature  to  use  a  subscript  g  for  the 
Taylor  microscale,  rather  than  a  subscript  z;  however,  due 
to  word  processing  limitations,  a  subscript  z  had  to  be 
used.  Xz  is  a  measure  of  the  average  dimension  of  the 
smallest  eddies  in  the  flow  field  [11,  p.  41].  In  order  to 
determine  the  size  of  Xz  in  terms  of  the  channel  width  D,  it 
is  necessary  to  introduce  the  equation  for  the  time  averaged 
turbulence  kinetic  energy.  This  equation  can  be  derived, 
with  much  work,  by  multiplying  equation  (2.2)  by  v  and  time 
averaging,  and  then  subtracting  from  this  the  equation 
obtained  by  multiplying  the  time  average  of  equation  (2.2) 
by  the  time  averaged  velocity  <v>  (not  to  be  confused  with 
<v> n  from  equation  (2.15)).  See  Hinze  [11,  pp.  68-72]  for  a 
derivation  in  cartesian  coordinates.  The  result  of  this  is 
the  equation  for  the  change  in  kinetic  energy  of  turbulence 
per  unit  mass  per  unit  time: 

0 . 5d ( <q2 > ) /dt  =  -V*<v’(p’/p  +q2/2)>  +  V* [ 2v<s «v ' > ] 

-  2 <  s  : Vv ' >  -  <v ?  v ' > : VV  ( 3 .  10) 

where  <q2>=<u2>+<v2>+<w2>  and  the  instantaneous  velocity  is 
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v=V+v '  (note:  <v>=V  and  v'=(u,v,w));  s  is  the  fluctuating 
stress  tensor  s  =  (Vv'  +  (Vv')t)/2;  and  p'  is  the 
fluctuating  component  of  pressure  i.e.  p  =  P  +  p'  (note: 
<p>=P) .  Despite  their  formidable  appearance,  each  of  the 
terms  in  equation  (3.10)  can  be  associated  with  a  physical 
process  [11,  p.  72].  According  to  Hinze  [11,  p.  72],  the 
last  term  in  the  equation  is  the  turbulence  production  term, 
and  the  second  last  term  is  the  rate  of  loss  of  <q2>  to  heat 
due  to  viscous  dissipation.  For  a  flow  in  which  all  time 
averaged  quantities  are  constant  in  time  (stationary 
turbulence),  we  can  make  the  local  equilibrium  assumption 
that  the  production  of  turbulence  is  equal  to  the 
dissipation  of  turbulence  [11,  p.  75].  In  a  channel,  the 
velocity  gradient  9U/9y  is  dominant  so  the  production  can 
be  written  [  14  ]  : 

p  =  -<uv>9U/9y  (3.11) 

The  viscous  dissipation  is  largely  due  to  eddies  which 
have  size  Xz  or  smaller,  and  these  eddies  are  formed  from 
much  larger  eddies.  By  the  time  the  eddies  have  reached  the 
size  Xz  they  are  no  longer  dependent  on  the  mean  flow  con¬ 
ditions  and  so  are  essentially  locally  isotropic.  Thus, 
because  most  dissipation  occurs  at  scales  that  are  locally 
isotropic,  we  can  use  the  exact  expression  for  isotropic 
turbulence  for  the  dissipation  e  [11,  p.  219]: 

e  =  ^5v<u2>/^h2z  (3.12) 

By  defining  A  as  the  average  size  of  the  energy  containing 
eddies  (i.e.  the  largest  eddies  in  the  flow),  then  e  is 
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governed  by  the  supply  of  turbulence  kinetic  energy  <u2>/2 
and  the  time  scale  ^/<u2>1/2  of  the  large  eddies.  Taking  e 
as  the  ratio  of  these  two  quantities,  then  [11,  p.  225]: 

e  =  A<u2>3/2//  (3.13) 

where  A  is  a  constant  near  unity.  We  can  also  estimate  the 
production  p  using  a  boundary  layer  approximation  for 
equation  (3.11)  since  9U/3y  ~  U02/D,  and  u  *  v  ~  <u2>1/2: 

p  ~  <u2>U0/6  (3.14) 

where  6=D/2.  Using  our  local  equilibrium  assumption  we  can 
equate  p  from  equation  (3.14)  with  e  from  equat i on ( 3 . 1 3 ) ,  so 
that  we  obtain  (as  in  [14]): 

4/b  '  <u2>1/2/U0  (3.15) 

Equating  equation  (3.12)  to  (3.13)  we  also  have  that: 

i/\z  =  0.258A1/2(<u2>1/2-^/*>)  1/2  (3.16) 

By  using  equations  (3.15)  and  (3.16)  we  have  the  desired 
result  for  Xz  in  terms  of  known  quantities: 

(D/2)/Xz  =  0.23(UoD/2p) 1/2  (3.17) 

where  a  value  of  A=0.8  has  been  used  in  equation  (3.13)  (see 
[11,  p.  248,  p.  255]).  Note  that  equation  (3.17)  was 
derived  under  the  assumptions  of  local  equilibrium  of  pro¬ 
duction  and  dissipation  and  for  turbulence  which  is 
isotropic  at  the  scale  Xz. 

One  can  also  define  a  third  eddy  size,  r\,  which  is 
associated  with  the  size  of  the  eddies  which  are  on  the 
verge  of  being  completely  dissipated.  This  is  the  smallest 
eddy  size  and  is  called  the  Kolmogoroff  scale.  r\  satisfies 
[ 1 1 , p.225] : 


43 


X  2  /  77  =  (15)  1/4*(<u2>1/2Xz/^)  1/2 


3.4  Analysis  of  distance  to  reach  steady  state  velocity 

profile 

In  modelling  channel  flow  numerically,  the  initial 
velocity  condition  of  the  channel  was  chosen  for  simplicity 
to  be  a  linear  profile  i.e.  the  along  stream  velocity  was 
initially  set  to  be  a  linear  function  of  distance  from  the 
wall  throughout  the  entire  channel.  Before  attempting  to 
simulate  the  channel,  it  is  thus  necessary  to  develop  a 
rough  estimate  for  the  time  required  for  this  velocity 
profile  to  transform  into  a  curve  resembling  the  fully 
developed  profile.  For  simplicity,  a  constant  velocity 
profile  not  a  linear  one  is  assumed  in  the  following 
derivation.  Note  that  a  linear  velocity  profile  fit  to  the 
experimental  data  would  obviously  approach  the  experimental 
profile  quicker  than  that  predicted  by  the  following  con¬ 
stant  velocity  profile  derivation. 

Assume  in  the  computational  grid  cell  nearest  the  wall 
that  the  viscous  shear  stress  term  p92U/9y2  is  about  equal 
to  the  Reynolds  stress  term  in  the  x-momentum  time  averaged 
equation 

pdU/dt  =  -9P/9x  +  p92U/9x2  +  p92U/9y2  -  p9(<uv>)/9y 

(3.18) 

Note  that  9U/9t=0  so  that  it  can  be  included  in  the  equation 
as  part  of  the  Lagrangian  derivative  without  any  trouble,  as 
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long  as  dU/dt  is  always  interpreted  as 
dU/dt=U9U/9x  +  V9U/9y.  Using  a  boundary  layer 
approximation,  then  certainly  9 2U/9x 2 <<29 2U/9y 2 .  In  order 
to  determine  when  9P/9x<<2p9 2U/9y 2  let  us  say  that  U  in  the 
first  cell  satisfies  U  ~  U0/2  and  U0  is  about  the  same  as 
<U>S,  the  average  value  of  U  across  the  stream;  then  the 
condition  that  9P/9x<<2p9 2U/9y 2  is: 

Xp<U>  2 /4D  <<  2pU/A2 

where  equation  (3.6A)  has  been  used  to  obtain  9P/9x  and  a 
boundary  layer  approximation  has  been  used  for  92U/9y2  with 
the  width  of  a  cell  A  as  the  length  scale.  Since 
U  "  U0/2  “  <U>s/2,  then  this  can  be  written  as 
XA/D  «  2p/UpA 

By  introducing  the  dimensionless  cell  width  A+=(u*)A/v,  then 
equation  (3.5)  can  be  written: 

UpA/p  =  A" ( 2/X ) 1 / 2 

so  that  the  condition  for  9P/9x<<2p9 2U/9y 2  is: 

A+  «  [ (X/2) 1 '  2 ]D/A  (3.19) 

X  is  always  less  than  0.04  for  turbulent  flow  [1,  p.  613], 
and  D/A  ~  30  (i.e.  the  channel  is  30  grid  cells  wide),  so 
that  equation  (3.19)  will  be  satisfied  for  A+<<200  or  so, 
i.e.  the  first  cell  should  not  extend  much  further  than  the 
buffer  layer  so  that  A+<30.  With  these  assumptions  equation 
(3.18)  can  be  written  for  the  cell  next  to  the  wall  as: 

pdU/dt  =  2p9  2U/9y 2  (3.20) 

This  equation  is  valid  throughout  the  first  cell  and  so  can 
be  integrated  over  the  volume  of  the  cell.  The  right  side 


can  be  expressed  as  a  surface  integral  since 

JV2UdV  =  /V*(VU)dV  =  /(VU)*dS  so  that  we  can  write,  for  a 

uniform  grid  spacing  and  assuming  3U/9x<<9U/9y : 
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pi (dU/dt)dV  &  p(-9U/9y | y=0  +  9U/9y|y=y1)A  (3.21) 

where  yl  is  the  y  position  of  the  point  midway  between  the 
first  and  second  cells  from  the  wall.  If  we  now  approximate 
the  left  hand  side  of  (3.21)  by  assuming  that  dU/dt  is  con¬ 
stant  at  dU1y/dt  over  the  cell,  as  it  would  be  in  a  finite 
difference  approximation,  then  p/( dU/dt )dv  =  pA2d.U,/dt, 
where  A2  is  the  volume  of  a  cell  of  unit  depth.  If  we  also 
assume  that  U  has  the  value  (U0+Ui)/2  in  the  second  cell  and 
approximate  the  derivatives  of  the  right  hand  side  of 
equation  (3.21)  with  finite  differences,  then  we  have: 

pA2dU i /dt  =  m(-5U!  +  U0)  (3.22) 

Now  dU/dt=9U/9 t+U9U/9x+V9U/9y ;  but  V=0  for  a  channel,  and 
9U/9t=0  by  definition,  so:  dU/dt=U9U/9x .  Making  the  substi¬ 
tution  £=x/U,  then  dU/dt=9U/9£  and  the  solution  to  equation 
(3.22)  is: 

U,U)  =  U0(4exp(-5p£/pA2 )  +1  )/5  (3.23) 

The  decay  constant  T  for  Ut(£)  is  thus: 

T  =  pA2/5p  (3.24) 

and  Ui  will  have  reached  a  value  0.5Uo  when  £=T.  Thus,  the 
downstream  distance  the  fluid  must  travel  before  the  value 
of  the  velocity  in  the  finite  difference  cell  next  to  the 
wall  changes  from  its  initial  value  of  U0  to  the  value 
0.5Uo,  because  of  viscous  drag  from  the  wall,  is:  xt=0.5Uor 


l  .e . 
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x  t  =  UoPA2/10/u  (3.25) 

This  provides  us  with  a  rough  estimate  of  what  distance  we 
would  expect  a  constant  velocity  profile  to  take  to  develop 
into  a  reasonably  sheared  velocity  profile  and  what  this 
distance  xt  depends  on. 


CHAPTER  IV 


Numerical  Simulation  of  Channel  Flow 

The  full  compressible  fluid  equations  (2.2)  and  (2.6) 
written  in  two  dimensional  cartesian  coordinates  x  and  y  are 
listed  in  Appendix  A.  These  were  modelled  on  the  University 
of  Alberta's  Amdahl  5860  computer  using  a  fluid  particle 
approach . 

4.1  Boundary  and  initial  conditions 

The  channel  geometry  simulated  was  that  of  a  fluid 
entering  a  channel  with  a  linear  velocity  profile,  as 
mentioned  in  section  3.4.  Obviously,  transition  to  a  fully 
turbulent  velocity  profile  could  not  be  expected  in  the 
short  length  channels  which  computation  times  demanded.  A 
typical  ratio  of  length  to  width  for  a  computational  run  was 
four,  whereas  one  would  not  expect  steady  state  flow  for 
much  longer  lengths  [11,  p.  709].  However,  by  the  addition 
of  a  random  component  to  the  entrance  fluid  velocities, 
enhanced  transition  to  turbulence  can  be  obtained 

[11,  p.  600] . 

As  an  initial  condition,  the  entire  channel  was  given  a 
linear  velocity  profile  i.e.  the  alongstream  velocity  was  a 
linearly  increasing  function  of  distance  from  the  wall.  The 
particles  were  given  velocities  which  yielded  a  Gaussian 
velocity  distribution  about  the  desired  initial  linear 
velocity  profile.  This  provided  an  initial  flow  that  had  a 
degree  of  random  motion.  A  standard  deviation  of  <U>S/10 
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was  used  for  the  Gaussian  distribution,  where  <U>S  is  the 
average  value  of  U  across  the  channel.  The  fluid  entering 
the  channel  was  also  given  a  linear  velocity  profile  with  a 
small  random  component  added.  The  particles  were  all  given 
the  same  initial  temperature  T0  and  the  number  of  particles 
per  cell  was  a  constant  n0,  so  that  the  pressure  and  density 
was  initially  constant  everywhere.  The  particles  were 
placed  randomly  within  each  cell. 

Since  the  channel  is  symmetric  about  its  centerline, 
only  one  half  of  the  channel  was  modelled.  The  channel 
geometry  and  ghost  points  are  depicted  in  figure  4.1.  The  x 
and  y  axes  are  offset  from  their  conventional  position  for 
convenience  of  defining  which  cell  a  particle  is  in.  On  the 
line  y=D/2  (j=Ny+2)  the  ghost  points  were  given  the  values 

( x  ^  i  ,  j  ( x  )  i  ,  j  —  1  r  (  V  y  )  '  .  j  ““  ,  j  -  1  '  P  i  ,  j  P  i  ,  j  —  1  » 

Ti(j=Tiij_1  because  the  full  channel's  centerline  lies  at 
y  =  D/2  -  A/2.  On  the  lower  boundary  v  must  be  zero  to 
satisfy  the  usual  no  slip  condition  at  the  wall.  As  well, 
T=T0 ,  and  8p/3y=0  here,  so  that  the  ghost  points  on  the  line 
y=-A  must  satisfy  vi(j=-vi(j+1,  Ti(j=2To-Ti)j+1  and 
p i  (  j  =p i  ,  j  +  i f  where  j=1.  The  particles  are  all  specularly 
reflected  from  the  y=-A/2  and  y  =  D/2  -A/2  boundaries.  The 
ghost  points  on  the  line  x=-A  were  given  values  which 
yielded  Neumann  boundary  conditions  for  v,  T  and  p  i.e. 

Xi.j=Xi*i,j'  T i  ( j =T i  +  i , j  f  Pi,j=Pi  +  i,j/  where  i  =  1. 

Similarly,  Neumann  conditions  were  used  for  v  and  T  for 
ghost  points  on  x=L;  however,  for  pressure  at  x=L  it  was 
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Figure  4.1  Ghost  point  configuration  for  channel  simulation. 
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found  that  special  consideration  had  to  be  given,  as 
described  in  the  following  section. 

4.2  Open  boundary  considerations  for  channel  flow 

4.2.1  Streaming  instability 

At  the  downstream  boundary  of  the  computational  window 
a  growing  instability  in  the  flow  arose  unless  proper  care 
was  taken  in  applying  the  pressure  boundary  condition  here. 
This  instability  has  a  very  physical  explanation;  if  the 
average  number  of  particles  in  a  cell  is  n0,  then  after  a 
particle  crosses  the  downstream  boundary  of  the 
computational  window,  the  cell  i,j  which  the  particle  left 
will  have  a  lower  pressure  in  general,  since  it  has  only 
ni(j=n0-1  particles  per  cell.  Thus,  the  particles  whose 
x-positions  lie  in  the  last  two  cells  of  the  channel  will 
feel  a  positive  x-direction  acceleration  from  the  -9p/3x 
force.  The  particles  thus  reach  the  end  of  the  channel  with 
a  higher  velocity  than  previous  particles.  If  a  Neumann 
condition  is  applied  to  the  pressure  at  the  boundary  here, 
then  the  (-9p/3x)  assigned  to  particles  in  the  last  two 
cells  of  the  channel  varies  linearly  between  zero  and  the 
positive  value  of  (pi  -  •!  ,  j~pi  ,  j  )/A  (using  a  modified  Vp 
differencing  with  pi  -  2  ,  j  “Pi  - 1  ,  j  and  p  i  _ ,  ,  j  >p  ■,  ,  j  ) .  However, 
particles  are  now  leaving  the  right  hand  side  of  cell  i,j  at 
a  faster  rate  than  they  enter  the  left  hand  side,  since 
particles  are  accelerated  across  the  last  two  cells.  Thus 
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ni(j  will  drop  faster  than  nj_1(j,  producing  a  larger 
(-dp/3x)  so  that  the  particles  leave  the  channel  with  an 
even  higher  velocity.  One  can  easily  see  that  this  particle 
accelerating  phenomena  will  work  its  way  upstream,  resulting 
in  a  completely  nonphysical  growing  loss  of  particles  out 
the  downstream  boundary.  This  undesirable  situation  can  be 
remedied  by  "clipping"  the  pressure  in  the  last  few  cells  of 
the  window.  By  not  allowing  the  pressure  to  drop  below  a 
certain  specified  value  in  these  cells,  the  above 
instability  cannot  develop.  That  is,  if  the  density  and 
temperature  are  such  that  the  actual  pressure  of  the  cell  is 
less  than  the  clipping  value,  then  for  all  finite 
differencing  the  clipped  value  was  used  instead  of  the 
actual  pressure.  It  was  convenient  to  use  equation  (3.3)  to 
obtain  an  estimate  of  the  pressure  drop  expected  across  a 
length  of  channel.  A  pressure  3.5%  of  p0  less  than  that 
predicted  from  (3.3)  was  then  applied  as  the  clipping  value 
in  order  to  minimize  unnecessary  clipping  of  turbulent 
pressure  variations.  To  ensure  that  the  instability  didn’t 
develop,  clipping  was  done  in  the  second  last  cell  as  well, 
but  with  a  lower  minimum  pressure  so  that  the  mininum 
permissable  pressure  in  the  second  last  cell  was  set  at 
about  2%  of  p0  less  than  that  of  the  last  cell. 

4.2.2  Buffer  zones 

The  application  of  pressure  clipping  removes  the  growth 
of  the  above  instability,  however  there  can  still  be  a 
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nonphysical  effect  resulting  from  particles  crossing  the 
window  boundaries.  At  the  downstream  boundary  this  effect 
is  the  sudden  appearance  of  a  false  pressure  drop  in  the 
cell  the  particle  left  (as  long  as  p>pd  where  pd  is  the 
clipping  pressure).  If  at  the  upstream  boundary,  then  the 
pressure  will  rise  suddenly  in  a  cell  where  a  particle  is 
introduced  into  the  window.  Both  of  these  anomalies  are  a 
result  of  suddenly  ignoring  a  particle's  contribution  to  the 
densities  in  the  cells  near  the  boundary  once  a  particle  is 
beyond  the  window.  The  solution  to  this  problem  was  the 
addition  of  buffer  zones  at  each  end  of  the  channel. 
Particles  in  the  buffer  zones  were  left  to  coast  with  no 
fluid  forces  acting  on  them.  Particles  were  created  at  the 
upstream  boundary  of  the  upstream  buffer.  Once  a  particle 
crosses  the  downstream  side  of  this  buffer  it  enters  the 
window  carrying  the  entrance  velocity  and  temperature  con¬ 
ditions.  However,  at  the  same  time,  another  particle  is 
created  at  the  upstream  side  of  the  buffer  so  that  in  a 
short  time,  after  this  particle  crosses  the  width  of  the 
upstream  buffer,  it  too  will  enter  the  window.  The  random 
component  of  the  entrance  velocity,  mentioned  earlier,  was 
obtained  by  giving  the  upstream  buffer  particles  the 
Gaussian  velocity  distribution  (which  was  also  applied  as  an 
initial  condition  to  the  window  particle  velocities).  The 
upstream  buffer  particles  were  also  given  a  linear  velocity 
profile,  in  order  to  maintain  the  linear  entrance  velocity 
prof i le . 
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When  a  particle  crosses  the  downstream  boundary  it  is 
placed  in  the  downstream  buffer  where  it  coasts  until  it 
crosses  the  downstream  boundary  of  the  downstream  buffer, 
after  which  it  is  dropped  from  memory.  The  density  of  the 
cells  in  the  channel  are  then  dependent  on  the  value  of  ni(j 
at  grid  points  inside  the  buffers.  In  this  way  particles  do 
not  cause  artificial  shocks  as  they  exit  and  enter  the 
window  but  are  treated  in  the  same  manner  as  a  particle 
inside  the  window  crossing  an  interior  cell  interface. 

It  is  also  convenient  to  calculate  the  value  of  the 
density  in  the  cells  of  the  downstream  buffer  which  are 
adjacent  to  the  window.  One  can  then  calculate  the  pressure 
in  these  cells  and  use  these  values  as  the  boundary  con¬ 
dition  pressures  for  the  last  cells  of  the  window.  For 
simplicity,  the  value  of  <T>n,  from  the  temperature  equival¬ 
ent  of  equation  (2.18),  was  used  for  the  temperature  in  the 
buffer  cells  in  which  pressure  (p=pRT)  was  calculated, 
rather  than  calculating  T(n)  from  the  temperature  version  of 
equation  (2.15). 

Both  buffers  were  given  an  initially  constant  value 
(T0)  for  the  particle  temperatures  and  initially  contained 
n0  particles  per  cell.  Both  buffers  were  made  long  enough 
so  that  a  particle  at  the  ends  furthest  from  the  window 
boundary  of  either  buffer  would  contribute  an  amount  just 
within  machine  precision  to  the  density  of  cells  inside  the 
window  at  the  window-buffer  boundaries.  For  a  Gaussian  mass 
density  profile  of  width  a2=A2/V2,  this  required  buffers 


54 


about  four  cells  long  for  single  precision  on  the  Amdahl 
5860.  The  two  buffers  were  also  given  Gaussian  particle 
velocities  with  standard  deviation  1/10  the  initial  total 
average  fluid  velocity. 

4.3  Other  considerations 

The  necessity  of  using  the  Vp  differencing  explained  in 
Chapter  2  was  mentioned  there.  It  should  be  noted  that 
using  the  modified  Vp  differencing,  described  in  Chapter  2, 
yielded  values  of  <u2>1/2/U0  which  were  as  small  as  one 
tenth  the  values  obtained  using  conventional  Vp 
differencing.  This  is  a  significant  improvement  for  the  NGP 
approach,  putting  it  on  the  same  level,  in  terms  of  random 
particle  noise  ("thermal  heating"),  as  the  SUD  grid  particle 
assignment  and  area  weighting  methods  described  in  [7]. 

4.3.1  Gaussian  averaging  of  v  and  T 

It  was  found  that  in  order  to  best  obtain  a  smooth 
fluid  picture  from  the  particles,  a  one  dimensional 
averaging  in  the  x-direction  was  desirable.  Previous  NGP 
codes  [2]  have  used  an  arithmetic  averaging  over  three 
cells;  this,  however,  reduces  resolution  considerably.  To 
obtain  a  local  smoothing  without  a  significant  reduction  in 
the  effective  grid  size,  a  one  dimensional  density  weighted 
Gaussian  averaging  with  width  a2=A2//2  was  used  on  velocity 
and  temperature.  That  is,  the  values  of  velocity  and  tem¬ 
perature  from  equation  (2.15)  (and  it's  temperature 
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equivalent)  were  averaged  in  the  x-direction  with  a  Gaussian 
weighting  before  being  used  in  the  finite  differencing  to 
obtain  L  and  E  in  equations  (2.16)  and  (2.17).  The  form  of 
the  Gaussian  averaging  written  for  the  fluid  velocity  at 
jTi  ,  j  at  time  t=nAt  is: 

(v{n))i  j  =  ( £_  i  ,  j-r*,,  j  )  (v(  n  5  )m,  jp£(  j  ]/ 

[Im/mof  (r  i  ,  j-r*, ,  j  )p£  ,  j  ]  (4.1) 

where  f  is  defined  in  equation  (2.8),  and  the  sums  are  over 
all  grid  points  in  the  jth  row  except  those  which  have  no 
particles  in  them,  the  mO  indicating  these  empty  cells. 

Empty  cells  are  excluded  from  the  averaging  because  v<n)  and 
T(n)  cannot  be  defined  in  an  empty  cell  using  an  equation  of 
the  type  given  by  (2.15).  The  value  of  (v(n))m,j  is  that 
from  equation  (2.15),  and  the  value  of  averaged  velocity 
(v(n))|  j  from  (4.1)  is  then  used  in  place  of  (v(n))i,j  to 
obtain  the  finite  differencing  in  the  L  and  E  of  equations 
(2.13)  and  (2.14).  For  a2=AJ/V2  (a=0.84A),  this  averaging 
uses  a  factor  f(0)=1;  a  factor  f(A)=0.49  for  cells  adjacent 
to  cell  i,j;  and  a  factor  f(2A)=0.06  for  cells  two  cells 
away,  compared  to  f(0)=f(A)=1  for  a  horizontal  arithmetic 
averaging  over  adjacent  cells.  The  averaging  is  thus  not  as 
large  scale  as  the  arithmetic  case. 

Without  the  Gaussian  averaging,  a  few  empty  cells 
developed,  while  virtually  none  existed  with  the  averaging. 
In  addition,  the  averaging  helped  reduce  the  random  noise 
caused  by  the  NGP  approximation  mentioned  in  section  4.5. 

The  removal  of  the  Gaussian  averaging  for  velocity  and 
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temperature  produced  turbulent  intensities  about  three  times 
that  obtained  for  a2=A2/\Z2,  in  an  early  run  with  am=at=0.1 
(am  and  at  are  from  equations  [2.16]  and  [2.17];  note  this 
is  too  large  a  value  of  am  for  channel  flow  -  see  section 
4.5) 

It  should  be  mentioned  that  in  the  computational 
versions  of  equations  (4.1)  and  (2.7),  only  cells  for  which 
f  was  within  machine  precision  were  used  in  calculating  the 
sums,  since  cells  further  away  would  have  insignificant  con¬ 
tributions.  The  subroutines  for  the  calculation  of  density 
from  equation  (2.7),  and  temperature  and  velocity  from 
equations  of  the  form  of  (4.1),  were  done  on  a  number 
crunching  array  processor  (FPS  190L)  attached  to  the  Amdahl 
5860. 


4.3.2  Empty  cells 

With  the  modified  Vp  differencing  described  in  Chapter 
2,  very  few  empty  cells  were  produced  in  any  of  the  runs 
(using  conventional  Vp  differencing,  one  in  ten  cells  was 
empty;  while  with  the  modified  Vp  there  were  usually  no 
empty  cells).  However,  the  particle  approach  described  in 
this  work,  which  obtains  the  fluid  density  based  on  the 
particle  positions  on  a  uniform  grid,  would  not  lend  itself 
well  to  modelling  a  system  with  order  of  magnitude  density 
variations,  since  this  would  require  some  cells  to  have  far 
more  particles  than  necessary  for  good  accuracy,  adding 
unnecessary  expense  to  a  run;  while  some  cells  would  have 
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only  a  few  or  even  no  particles,  resulting  in  poor  resol¬ 
ution  of  the  fluid  here.  A  variable  grid  might  be 
introduced,  but  this  would  add  problems  of  its  own.  For 
those  interested  in  modelling  a  system  that  might  have  empty 
cells,  or  "holes"  as  we  have  called  them,  it  is  critical  to 
note  that  holes  cause  various  problems  in  the  conservation 
of  fluid  quantities  when  using  flux  conservative 
differencing.  This  is  because  there  are  no  particles  in  the 
hole  which  can  cancel  the  flux  conservative  differencing 
from  neighboring  cells.  Consider,  for  example,  the 
( d 2T/9x 2 ) i , j = (T j + t , j +T j _ i  j -2T j  , j ) /A2  term  appearing  in  E 
from  equation  (2.17).  (Actually  the  temperature  used  here 
is  that  from  the  temperature  equivalent  of  (4.1)  but  the 
circumflexes  will  be  left  off  for  convenience.)  With  no 
empty  cells,  then  ( 3  2T/3x  2 ) i+1 ( j=(Ti+2, j+Tj , j-2Ti+1 ( j  ) /A2 
and  (  3  2T/9x 2 )  i  _  !  ;  j  =  (T  j  ,  j +T  i  .  2  ,  j  ~2T  j  _  ,  j  )/A2  ,  so  that  the  sum 
on  the  right  hand  side  of  equation  (2.21)  would  properly 
vanish.  However,  if  cell  i,j  is  empty,  then  there  are  no 
particles  which  will  yield  a  value  of  (32T/3x2)iij  in  the 
inner  sum  of  the  temperature  equivalent  of  equation  (2.19). 
Thus,  the  -2Ti(j  term  in  (32T/3x2)i(j  (note  that  with 
Gaussian  averaging,  the  fluid  temperature,  velocity  and 
density  are  defined  even  in  an  empty  cell)  will  not  cancel 
the  Ti,j  term  of  ( 3 2T/3x 2 )  j  +  ,  f  j  and  the  T^j  term  of 
(3 2T/3x 2 ) j _ i , j ,  so  that  the  right  hand  side  of  equation 
(2.21)  is  nonzero,  and  we  no  longer  have  energy 
conservation.  This  situation  obviously  occurs  for  all  the 
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terms  in  both  the  momentum  and  energy  equations  which  are 
differenced  flux  conservatively,  and  which  should  thus 
vanish  when  summed  over  the  grid. 

For  all  terms  except  the  Vp  term  of  the  momentum 
equation  this  situation  can  be  partially  remedied  by 
applying  a  Neumann  condition  for  the  values  of  temperature, 
velocity,  and  pressure  at  the  interface  of  the  empty  cell. 
For  example,  instead  of  using  the  value  of  Ti(j  in  the 
expression  for  ( 9 2T/9x 2 ) i + i ,  j ,  one  would  use  the  value  of 
Ti+1)j,  and  for  ( 9 2T/9x 2 ) f _ , ( i  one  would  use  Tj_1(j  in  place 
of  T s  j .  In  this  way  all  92A/9x2,  92A/9y2  terms  will  vanish 
exactly  when  summed  over  the  grid.  Unfortunately,  the 
single  derivative  terms  and  the  mixed  derivatives  will  not 
exactly  vanish  using  this  scheme,  since  their  flux 
conservative  finite  difference  forms  are  based  on 
cancellation  over  cells  which  are  two  grid  points  apart. 

For  example,  (  9T/9x )  •,  ,  }  =  (T  j  +  ,  ,  j  -T  j  .  ,  f  j  )/2A  is  normally 
cancelled  by  a  -Ti+1(j  from  ( 9T/9x ) , . 2 , j  and  a  Tj_1(j  from 
( 9T/9x ) i . 2 ,  j .  If  cell  i  —  2 , j  is  empty  then  there  will  be  no 
T  j  _  i  (  j  from  (  9T/9x )  j  _  2  ,  j  ,  but  there  will  be  an  extra  -T  j  _  •,  ,  i 
by  applying  a  Neumann  condition  for  (9T/9x) i . i ,  j . 

Similarly,  we  will  accumulate  2T  j  _  3  j  on  the  other  side  of 
the  hole.  If  2T j  _  3  f  ]  =2T j  _  •,  (  }  ,  then  the  sum  will  vanish,  but 
in  general  this  will  not  be  exactly  true.  This  difficulty 
is  illustrated  schematically  in  figure  4.2.  In  the  top  left 
figure,  a  flux  conservative  finite  differencing  of  92A/9x2 
will  vanish  exactly  when  summed  over  the  grid;  similarly  for 
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9A/9x  as  in  the  top  right  figure.  However,  when  a  hole 
appears  in  cell  i,j  (as  depicted  by  the  shaded  region),  then 
a  term  -A j . , ( j +2A j ( j -A j + , # j  (as  depicted  in  the  middle 
figure),  will  be  present  in  a  sum  of  92A/9x2  over  the  grid, 
so  that  the  differencing  is  no  longer  conservative.  By 
applying  a  Neumann  condition  at  the  empty  cell  interface, 
the  problem  is  corrected,  as  in  the  bottom  left  picture. 
Unfortunately,  9A/9x  remains  nonconservative.  In  the  middle 
right  figure,  we  see  an  extra  A  j  _  ,  t  j -A  ■,  +  ,  i  j  appears  with  no 
Neumann  condition;  while  applying  the  Neumann  condition 
gives  an  extra  2A  ■,  _  •,  ,  } -2A  j  +  1  ,  j  as  in  the  bottom  right 
figure.  It  is  clear  that  it  would  not  be  wise  to  use  a 
Neumann  condition  for  evaluating  Vp  in  cells  adjacent  to  a 
hole,  since  there  would  then  be  no  tendency  for  the  hole  to 
be  filled  in.  Note  that  over  a  long  time  period  with  steady 
state  conditions  one  would  expect  the  nonzero  terms  in  the 
sums  of  L  and  E  over  the  grid,  which  are  caused  by  empty 
cells,  to  average  to  zero. 

4.4  Constraints  on  variables 

The  quantities  involved  in  the  description  of  channel 
flow  must  satisfy  various  conditions  before  we  would  expect 
a  simulation  to  converge  to  the  correct  solution.  As 
mentioned  in  Chapter  1,  the  Courant  (CFL)  stability 
criterion  must  be  satisfied,  i.e. 

{ | v j  | mx  +  n/p&  +  ( C s  2  +v t  2 ) 1 7  2 } At/A  <  1  (4.2) 

where  |vj|mx  is  the  maximum  particle  velocity;  Cs  is  the 
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Figure  4.2  Cells  with  no  particles  in  them  introduce  a 
nonconservative  effect  into  the  normally 
conservative  flux  conservative  finite 
differencing.  A  plus/minus  sign  above  a  cell 
signifies  that  the  value  of  A  from  that  cell  is 
added/subtracted  to/from  a  sum  of  the 
differencing  over  the  grid.  (1)  No  empty  cells, 
(2)  empty  cells  -  no  corrective  action  in  the 
differencing  and  (3)  Neumann  treatment  of  empty 
cells. See  p.  59  for  a  further  explanation. 
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adiabatic  sound  speed  (for  an  ideal  gas  Cs=[7P/p]1/2  where  7 
is  the  ratio  of  specific  heats);  vt  is  the  Alfven  wave  speed 
(vt  2  =B 2 /p0  P  for  a  plasma,  where  B  is  the  magnetic  field 
strength  and  p0  is  the  permeability  of  free  space);  and  p/pA 
is  the  speed  at  which  the  velocity  is  diffused  by  viscous 
effects.  For  the  case  of  the  simulations  presented  here, 
p/pA<<Cs;  and  also  the  magnetic  field  fluctuations  are  not 
considered  in  the  equations  (i.e.  there  is  no  magnetic  field 
equation),  therefore  it  is  not  appropriate  to  include  vt  in 
equation  (4.2),  so  (4.2)  can  be  written  as: 

At  <  A/( | Vi  | m x  +C s )  (4.3) 

The  grid  size  A  is  also  limited  by  the  requirement  that  A  be 
less  than  or  of  the  same  size  as  the  eddy  size  Xz  from 
equation  (3.17),  i.e.  A/X2  <  1 ,  so  that  from  equation 
(3.17)  with  6=NyA=D/2,  where  Ny  is  the  number  of  cells 
across  the  half  channel,  then: 

0.23(NyUopA/p) 1 /2/Ny  <  1  (4.4) 

Since  we  are  interested  in  modeling  turbulent  fluid  flow,  it 
is  necessary  to  maintain  the  Reynolds  number  above  the 
turbulent-laminar  transition  region.  For  Re d =<U> sD/^<2300 , 
where  <U> s  is  the  mean  velocity  of  the  channel,  the  flow  is 
nearly  always  laminar  [1,  p.39].  The  change  to  turbulent 
flow  occurs  over  a  range  of  Reynolds  numbers  slightly  above 
this.  With  very  undisturbed  inlet  conditions,  the 
transition  to  turbulence  may  not  occur  until  Red>10,000; 
however,  it  was  mentioned  earlier  that  well  disturbed  inlet 
conditions  were  used.  A  Reynolds  number  of  Red=5333  was 
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thus  regarded  as  a  safe  lower  limit,  i.e. 

<U>SD /v  >  5333  (4.5) 

In  order  that  the  fluid  develop  a  velocity  profile  that 
isn’t  too  far  from  the  steady  state  profile  by  the  time  it 
reaches  the  end  of  the  window,  we  require  that  the  window  be 
at  least  about  twice  the  distance  xt  given  in  equation 
(3.25).  With  the  length  of  the  channel  given  by  L=NXA  then 
this  requirement  is  (letting  U0=<U>S): 

Nx A  >  <U>  s  pA2 /5m  (4.6) 

The  time  required  for  the  fluid  to  cross  the  window  and 
develop  the  correct  velocity  profile  must  also  be  within  a 
reasonable  number  of  time  steps  for  expense  reasons.  If 
this  number  of  time  steps  is  NTX ,  and  At  has  been  chosen  as 
a  constant  At=rA/(Cs+<U>s ) ,  where  |r|<1;  then  (3.25)  gives: 

NTX  =  ( C  s  +<U>  s  )  pA/5jur  (4.7) 

Since  the  data  of  Nikuradse  and  Donch  (figure  3.2)  is  for 
incompressible  flow,  then  the  velocity  of  the  flow  must  be 
low  enough  that  compressible  effects  do  not  become  signifi¬ 
cant.  From  [1,  p.  10]  this  requirement  is: 

U0  <  0 . 3C $  (4.8) 

where  U0  is  the  midstream  velocity.  Also,  the  length  of  the 
channel  should  be  such  that  it  is  at  least  about  ten  large 
scale  eddy  lengths  long  in  order  to  have  reasonable 
turbulent  modelling.  From  equation  (3.15)  with 
<u2> 1 ; 2/U0  ~0 . 1  as  mentioned  in  Chapter  3,  then  this 
requires : 

L/6  >  1 
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i.e.  the  length  to  width  ratio  of  the  window  should  be 
greater  than  1  or  so. 

To  satisfy  the  requirements  of  equations  (4.3)  to 
(4.8),  the  following  values  of  constants  were  chosen: 
r=0.2,  Cs=374  m/s,  Ny=16,  Nx=100,  p= 1.0  kg/m3, 

<U>  s  =  112.2  m/s,  m=1 .74169X10'“  kg/ms,  A=2 . 5875X1 0 ' “  m. 
This  yields  a  Reynolds  number  of  Re d  =<U> sD/i^  =  5333 .  <U>S  is 
constant  in  time  because  the  left  buffer  maintains  a  con¬ 
stant  rate  of  flow  into  the  computational  window.  One 
should  realize  that  at  this  low  Reynolds  number  the 
turbulence  is  probably  not  isotropic  at  the  scale  X2,  since 
the  large  scale  eddy  size  is  about  the  size  of  the  Taylor 
microscale,  based  on  equation  (3.16).  As  a  result,  eddies 
of  size  X2  are  likely  still  affected  by  the  eddies  of  size 
4,  which  are  not  isotropic  since  the  large  scale  eddies  are 
dependent  on  the  mean  flow.  Equation  (3.17)  is  thus  only  a 
very  rough  approximation  for  X2.  In  addition,  the  following 
other  parameters  were  chosen  as: 

p| t =o=100,000  Pa,  T0=348  K,  Cv=732  J/kgK, 
k=2 . 94X1 0 " 2  J/smK 

Also,  the  linear  velocity  profile  was  characterized  by  a 
velocity  of  81  m/s  at  the  wall  and  155  m/s  on  the  channel 
centerline.  Note  that  with  these  values  the  dimensionless 
cell  width  A*  ,  from  Chapter  3,  is  A+=11;  thus,  the  first 
cell  extends  only  into  the  buffer  layer  so  that  equation 
(3.19)  is  satisfied,  and  as  a  result  equation  (3.25)  is 
valid.  It  was  found  that  using  n0=m2  particles  per  cell, 
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where  m  is  an  integer,  was  desirable  for  random  placement  of 

the  particles.  In  this  way  each  particle  could  be  placed 

randomly  in  an  area  occupying  1/m2  of  a  cell.  Values  of  m, 

« 

as  well  as  am  and  a,  (from  equations  [2.16]  and  [2.17])  are 
discussed  in  the  next  two  sections. 

4.5  The  artificial  drag  term  and  noise  caused  by  a  NGP 

approximation 

It  was  initially  thought  that  the  values  of  am  and  at 
used  in  a  simulation  (am  and  at  are  the  artificial  drag 
coefficients  in  equations  [2.16]  and  [2.17])  would  not 
affect  the  results  of  the  simulation  in  any  significant  way, 
since  they  contribute  no  net  force  to  any  cell.  Of  course, 
if  too  small  values  are  used  then  there  is  a  large  amount  of 
mult i streaming  and  the  simulation  becomes  too  noisy  to  yield 
any  realistic  results.  However,  for  the  channel  simulation, 
it  was  found  necessary  to  have  am  such  that  the  intracell 
artificial  drag  forces  were  of  the  same  size  as  the 
intercell  molecular  viscous  forces.  It  is  not  known  exactly 
why  this  is  necessary,  but  it  is  reasonable,  since,  in  a 
sense,  the  artificial  drag  force  models  intracellular 
viscosity  so  that  intracell  viscosity  (caused  by  the 
artificial  drag  term)  should  be  the  same  as  intercell 
viscosity  (caused  by  the  actual  viscous  terms  appearing  in 
the  fluid  equations).  In  a  turbulent  boundary  layer  flow 
where  molecular  viscous  effects  are  expected  to  extend  to 
subgrid  length  scales,  using  a  multistreaming  term  that  is 
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much  larger  than  the  dominant  viscous  term  is  not  realistic 
and  causes  the  effect  of  intercell  viscous  forces  to  be 
swamped  by  the  effect  of  intracell  drag  forces,  resulting  in 
an  unrealistic  flow.  However,  in  a  laminar  flow,  large 
values  of  am  could  probably  be  used  without  any  difficulty 
since,  with  reasonable  grid  resolution,  there  should  be  no 
subgrid  velocity  variations,  so  that  the  particle  velocities 
should  be  very  closely  tied  to  the  fluid  velocity. 

If  too  large  a  value  for  am  was  used  in  the  turbulent 
channel  simulation,  various  incorrect  velocity  profiles 
resulted.  Although  the  temperature  equation  was  used  for 
stability  in  the  channel  code,  the  channel  flow  was  not  a 
thermal  boundary  layer  so  that  the  effect  of  large  at  on  the 
flow  was  not  destructive.  However,  in  a  flow  where  subgrid 
conduction  effects  exist,  one  would  probably  have  to  use  a 
value  of  at  which  yielded  intracell  (artificial)  conduction 
terms  which  were  the  same  size  as  intercell  (actual) 
conduction  terms.  In  the  gas  discharge  simulation,  the  tem¬ 
perature  profile  was  almost  completely  determined  by  Joule 
heating  of  the  fluid,  so  that  using  large  values  of  at 
(at=0.5)  was  also  not  destructive.  As  well,  for  the 
discharge  geometry,  large  values  of  am  were  acceptable, 
since  there  was  no  boundary  layer  flow  in  the  r,z  velocities 
and  the  6  velocity  was  essentially  laminar  because  of  the 
two  dimensional  nature  of  the  code. 

Let  us  now  develop  an  estimate  of  the  order  of  magni¬ 
tude  of  actual  viscous  particle  accelerations  and  of 
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particle  accelerations  caused  by  the  mult i streaming  term. 

The  actual  change  in  along  stream  velocity  that  a 
particle  undergoes  in  one  time  step  due  to  molecular  viscous 
effects  is:  ( p9 2vx/3y 2 )At/p ,  which  can  be  estimated  as: 
p<U> s At/pA2 .  The  change  in  a  particle's  velocity  due  to 
artificial  drag  can  be  estimated  as  follows.  If  we  say 
that,  at  time  nAt ,  a  particle's  velocity,  v|)"1/2,  differs 
from  <v>"  of  equation  (2.18)  by  <U>S 

i.e.  (yk'1/2  -  <v>")  ~  <U>S,  then  the  artificial  drag  term 
will  change  that  part icle ' s  veloc ity  by  am<U>s.  If  we  were 
to  have  the  effects  of  the  molecular  viscous  and  artificial 
drag  terms  about  equal  per  time  step,  then  this  would 
require  that  am<U>,  =  ju<U>sAt/pA2  so  that: 

am  -  pAt/pA2  (4.9) 

Note  that  the  effect  of  the  artificial  drag  term  on  a  par¬ 
ticular  particle  decreases  geometrically  in  time,  since  the 
magnitude  of  the  change  in  that  particle's  velocity  caused 
by  this  term  depends  on  |v£'1/2  -  vn  |  ,  which  is  decreasing 
in  time  because  of  the  effect  of  the  artificial  drag  term. 
Thus,  it  might  be  better  to  compare  the  two  terms  based  on 
how  many  time  steps  they  take  to  change  a  particle's 
velocity  by,  say,  <U>s/2.  For  the  molecular  viscous  term 
let  this  number  of  time  steps  be  K,  so  that 
K ( p<U> t At/pA2 )  =  <U>s/2  i.e.  K  =  pA2/2pAt,  or,  using 
At  =  rA/(Cs+<U>s)  as  in  obtaining  equation  (4.7),  then: 

K  =  p(Cs+<U>s )A/2pr  (4.10) 

(note  the  similarity  to  4.7).  For  the  artificial  drag  term, 
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let  us  say  that  a  particle's  velocity  is  initially  <U> s 
above  <v>n;  then,  the  number  of  time  steps,  M,  needed  to 
reduce  a  particle's  velocity  by  <U>i/2  is  found  from 
(1-am)M  =1/2,  so  that: 

M  =  -/n2//n(  1-0  (4.11) 

In  order  to  have  K=M,  then  we  must  have  am  such  that 

am  =  1  -  2 ‘ 1 (4.12) 
where  K  is  given  by  (4.10).  Equation  (4.9)  would  suggest  a 
value  of  am=0.0003  using  the  parameters  listed  in  section 
4.4,  while  equation  (4.12)  would  suggest  a  value  of 
am=0.0004.  If  values  of  am  much  larger  than  this  were  used, 
incorrect  velocity  profiles  resulted,  as  mentioned  earlier. 
For  example,  it  was  found  that,  although  noise  levels  were 
neglible  with  a  value  of  am=at=0.1  with  9  particles  per 
cell,  an  incorrect  velocity  profile  resulted,  as  shown  in 
figure  (4.3).  A  value  of  am=0.01,  n0=9  also  yielded  an 
incorrect  velocity  profile.  A  value  of  am=0.0015  was  used 
successfully  in  the  channel.  Unfortunately,  when  values  of 
am  this  small  were  used,  the  particle  velocities  took  on 
random  components  that  were  as  large  as  actual  turbulent 
fluctuations.  These  random  velocities  resulted  largely 
because  of  random  variations  in  pressure,  which  were  in  turn 
caused  by  random  variations  in  density  and  temperature.  The 
random  fluctuations  in  density  and  temperature  were  a  result 
of  the  discrete  nature  of  a  NGP  fluid  particle  method. 

Using  a  value  of  at=0.5  kept  temperature  fluctuations  to  low 
levels  and  also  helped  reduce  the  particle  velocity  noise 
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Figure  4.3  A  velocity  profile  that  results  for  channel  flow 

if  too  large  a  value  of  the  artificial  drag 
coefficient  am  is  used.  The  code  results  are 
indicated  by  the  triangles,  while  the  experi¬ 
mental  results  of  Nikuradse  and  Donch 
[12,  p.  482]  are  indicated  by  X  and  0  respect¬ 
ively. 
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levels,  however  a  way  of  reducing  the  random  density 
fluctuations  needs  to  be  developed;  that  is,  if  using  more 
particles  per  cell  is  not  an  acceptable  means  of  reducing 
noise.  These  random  density  variations  result  largely  from 
particles  flowing  across  cell  boundaries.  As  the  particles 
flow  across  the  grid,  the  number  of  particles  per  cell  ni(j 
fluctuates  randomly  because  there  are  a  finite  number  of 
particles  randomly  placed  in  a  cell.  Using  a  NGP  grid 
particle  assignment,  this  results  in  random  density 
variations  which  in  turn,  through  the  Vp  force,  yields  noise 
in  the  particle  velocities. 

To  quantify  this  effect  of  the  NGP  approximation  on 
noise  levels,  consider  a  simple  one  dimensional  case  where 
all  cells  have  m2  particles  in  them  (m  is  an  integer), 
except  cell  i+1,j,  which  has  ni+1(j  =  m2-6  particles  in  it. 
Here,  the  fluctuation  in  n  is  due  simply  to  a  finite  number 
of  randomly  placed  particles  flowing  across  a  grid;  it  is 
not  a  fluctuation  due  to  actual  fluid  turbulence.  With  a 
value  of  a2=A/i/2,  then  f=0.5  in  adjacent  cells  and  f=0  for 
other  cells,  so  that,  using  equation  (2.7): 

Pi  j  =  p o ( 1  -  6/4m2 ) ,  pi  +  i,j  =  p0(1  “  26/4m2 ) .  Thus, 
between  cells  i+1,j  and  i,j  a  Vp  force 

Vp  =  p06/4m2A  (4.13) 

develops  (ignoring  the  effects  of  a  modified  Vp 
differencing)  due  to  the  fluctuation  in  n.  In  one  time 
step,  the  change  in  a  particles’s  velocity  due  to  this  false 
Vp  force  is: 
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Av  =  VpAt/po  (4.14) 

If  the  velocity  of  the  fluid  in  cell  i,j  is  U  and  if  it 
takes  x  steps  for  a  particle  flowing  at  speed  U  to  cross  a 
cell  (x=A/UAt),  then  it  will  take  about  x6/m2  time  steps 
before  this  density  gradient  is  corrected  by  other  particles 
flowing  into  the  cell  i+1,j  from  cell  i,j.  This  is  because 
there  are  6  particles  in  a  6/m2  area  fraction  of  the  cell, 
so  that  6  particles  should  flow  into  cell  i+1,j  in  x6/m2 
time  steps  (assuming  square  cells).  In  two  dimensions, 
there  are  three  other  cells  adjacent  to  cell  i+1,j  from 
which  particles  could  also  flow  into  cell  i+1,j,  but  let  us 
assume  that  most  of  the  particles  move  into  cell  i+1,j  from 
the  upstream  cell,  i , j .  Using  equations  (4.13)  and  (4.14), 
then  in  x6/m2  steps,  the  change  in  a  particle’s  velocity 
would  be: 

x6Av/m2  =  p062/4p0Um4 

which  creates  a  turbulent  intensity  of  order: 

u ' /<U>  s  =  p062/4p0Um4  <U> ,  (4.15) 

where  u'=x6Av/m2.  This  suggests  that  the  noise  levels 
caused  by  particles  crossing  a  boundary  should  be  much  less 
for  runs  with  higher  values  of  m2 .  However,  the  above 
analysis  fails  to  account  for  the  fact  that,  for  m2  not  too 
large,  the  value  of  6  would  likely  increase  as  m2  increases 
(for  large  m2  this  would  not  be  true  since  one  would  be 
approaching  a  continuous  system),  because  the  fluctuation 
across  the  cell  boundaries  is  a  result  of  the  random 
placement  of  particles.  With  more  particles  randomly  placed 
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in  a  cell,  6  would  be  larger,  although  the  percentage 
fluctuation  in  the  number  of  particles  per  cell  would  be 
less.  If  we  say  that  6  *  m,  (4.15)  becomes: 

u ’ / <U>  s  *  p0/4p0Um2  <U>  s  (4.16) 

Using  the  values  of  p0 ,  Po,  <U> s  listed  earlier  and  assuming 
U=<U>S,  then  this  yields  noise  levels  of  u' /<\J> ,  =  0.23, 
0.12,  0.08  for  m2  =  9,  16,  25  respectively.  Equation  (3.9A) 
yields  an  analytic  turbulence  intensity  of  <u2 > 1 ; 2/<U> s=0 . 09 
so  that  the  noise  levels  are  of  the  same  order  as  the 
turbulent  intensities.  Computationally,  it  was  found  that 
approximate  values  of  <u2>1/2/<U>s  with  at=am=0.0015  were: 
0.35,  0.26,  0.25  for  n0=9,  16,  25;  subtracting  the  turbulent 
intensity  of  0.09  from  these  values,  then  noise  levels  of 
0.26,0.17  and  0.16  were  observed  for  n0=9,  16,  25  respect¬ 
ively.  Comparing  the  values  of  0.26,  0.17  and  0.16  from  the 
code  to  the  values  of  0.23,  0.12,  0.08  from  equation  (4.16), 
it  appears  that  the  noise  level  due  to  random  density 
fluctuations  does  not  decrease  quite  as  rapidly  as  1/m2  with 
increasing  m2  as  (4.16)  predicts,  although  equation  (4.16) 
is  only  an  order  of  magnitude  estimate.  Note  also  that  we 
have  not  accounted  for  the  noise  caused  by  temperature 
variations  mentioned  earlier.  Since  temperature  variations 
appear  directly  as  pressure  variations,  one  might  expect  the 
noise  caused  by  such  variations  to  be  of  the  same  order  as 
the  density  fluctuation  noise  levels.  Increasing  at  to  0.5 
removed  most  of  the  temperature  fluctuations  and  thus 
yielded  less  particle  velocity  noise.  With  at=0.5  and 
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am=0.0015,  runs  with  n0=9  and  25  gave  values  of 

<u2> 1 7  2/<Us>=0 .22  and  0.16  respectively,  considerably  lower 

than  the  values  of  0.35  and  0.25  found  with  at=0.0015. 

Because  the  particle  noise  is  largely  random,  its 
effect  on  time  averages  like  U  and  P  can  average  out  so  that 
good  results  can  be  obtained  for  these  quantities,  even  with 
low  values  of  m.  Physically,  however,  if  too  much  random 
mixing  occurs  across  the  channel,  then  the  velocity  profile 
is  diffused  into  a  constant  value  across  the  channel.  This 
was  in  fact  observed  in  earlier  runs  which  used  the  values 
am=a t =0 . 00 1 5 .  For  quantities  like  <u2>  and  <uv>,  noisy 
results  are  obtained  since  the  particle  noise  is  inseparable 
from  turbulence  in  these  averages  and  so  the  noise  does  not 
average  out. 

Although  the  particle  noise  caused  by  a  NGP 
approximation  interferes  with  obtaining  values  for  turbulent 
correlations  like  <u2>  and  <uv> ,  the  noise  does  not  have  a 
large  effect  on  the  fluid  dynamics,  as  long  as  noise  levels 
are  reasonable.  This  is  clear  when  one  considers  that  with 
at=0.5,  the  temperature  fluctuations  are  small  so  that  the 
major  source  of  the  noise  is  a  random  difference  in  the 
number  of  particles  in  a  cell.  These  random  differences 
last  only  a  fraction  of  the  time  required  for  a  particle  to 
cross  a  cell,  as  mentioned  in  deriving  equation  (4.15).  As 
well,  these  random  differences  are  likely  overcompensated 
for,  because  of  the  lag  time  between  the  random 
accelerations  and  the  random  velocities.  That  is,  if 
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particles  in  an  adjacent  cell  are  randomly  accelerated 
toward  a  cell  with  fewer  particles,  they  will  maintain  a 
higher  velocity  in  that  direction  even  after  the 
acceleration  stops,  resulting  in  too  many  particles  entering 
the  cell;  the  particles  are  then  accelerated  in  the  opposite 
direction  -  the  result  is  an  oscillation  of  particles 
between  adjacent  cells,  which  does  not  contribute  much  to 
the  fluid  motion. 

4.6  Results 

With  these  considerations  made,  good  results  were 
obtained.  A  comparison  of  the  time  averaged  along  stream 
velocity  with  experimental  data  is  shown  in  figure  4.4  for 
n0=25,  am=0.0015,  at= 0.5.  The  data  of  Nikuradse  and  Donch 
as  it  appears  in  figure  3.2  was  unnormalized  as  described  in 
section  3.1.  The  code  results  were  obtained  from  a  time 
average  over  the  last  100  time  steps  and  a  lengthwise 
spatial  average  over  the  section  x=C.8L  to  x=0.96L.  With 
this  averaging  not  all  turbulence  was  removed  from  U,  so 
that  the  code  results  oscillated  by  small  amounts  about  the 
experimental  data.  This  is  seen  in  figure  4.4  near  2y/D=0.3 
where  the  code  results  are  a  bit  lower  than  the  experimental 
data  and  at  2y/D=0.9  where  the  code  results  are  a  bit  higher 
than  the  experimental  data.  These  oscillations  occurred  at 
different  values  of  2y/D  for  different  time  steps.  With  9 
particles  per  cell,  the  velocity  profile  agrees  well  with 
experimental  results  by  600  time  steps.  Runs  up  to  1800 


Figure  4.4  Comparison  of  code  velocity  profile  (triangle 
with  experimental  data  of  Nikuradse  (x)  and 
Donch  (0)  from  [12,  p.  428]  after  4200  time 
steps  (4.47  ms ) . 
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time  steps  were  done  with  9  particles  per  cell;  the  code 
continued  to  maintain  profiles  as  in  figure  4.4  for  the  full 
length  of  these  runs.  With  25  particles  per  cell,  random 
noise  levels  were  lower  so  that  turbulence  required  a  longer 
time  to  develop,  since  turbulence  development  is  enhanced  by 
the  random  noise  (as  mentioned  in  section  4.1).  As  a 
result,  runs  with  25  particles  per  cell  required  a  longer 
time  to  develop  a  steady  state  velocity  profile.  For  n0=25, 
steady  state  fully  developed  velocity  profiles  were  reached 
by  3600  time  steps  and  were  maintained  for  the  remainder  of 
the  full  run  to  4200  time  steps.  The  initial  condition  and 
channel  entrance  profile  was  a  straight  line  profile  through 
the  points  2y/D=0.0,  U=81  m/s  and  2y/D= 1 . 0 ,  U=155  m/s,  as 
mentioned  earlier. 

A  plot  of  the  time  averaged  pressure  versus  along 
stream  distance  is  shown  in  figure  4.5  with  the  analytic 
pressure  profile  from  equation  (3.6A)  also  plotted.  The 
code  results  were  time  averaged  over  the  last  100  time  steps 
and  spatially  averaged  across  the  width  of  the  channel.  One 
would  not  expect  the  pressure  profile  to  be  correct  at  the 
entrance  to  the  channel  since  fully  developed  conditions  do 
not  exist  there. 

A  typical  snapshot  of  the  instantaneous  flow  in  the 
channel  is  shown  in  figure  4.6,  where  the  fluid  velocities 
at  all  grid  points  in  the  section  x=0.4L  to  x=0.7L  are 


shown . 
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Figure  4.5  Comparison  of  pressure  profile  (triangles)  to 
analytic  results  (straight  line). 


2Y/D 
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Figure  4.6  A  typical  instantaneous  flow  pattern.  This 

snapshot  is  for  3100  time  steps  (3.3  ms).  The 
arrows  are  scaled  to  the  maximum  velocity,  which 
is  186.6  m/s  (Note:  <U>,=112.2  m/s) 
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In  the  simulation  of  the  channel,  the  maximum  number  of 
25  particles  per  cell  was  limited  to  25  because  of  excessive 
computing  costs  required  if  more  particles  were  used.  With 
these  values  of  n0,  the  random  particle  noise  mentioned 
earlier  was  of  the  same  order  as  expected  levels  of 
turbulent  fluctuations,  so  that  turbulent  correlations  like 
<u2>  and  <uv>  did  not  agree  as  well  with  experimental  and 
analytic  results.  Code  results  for  the  shear  stress,  r,  as 
a  function  of  cross  stream  distance  were  usually  scattered 
about  the  straight  line  profile  expected  (equation  [3.7]). 
The  averaging  used  to  obtain  r  was  the  same  as  that  used  in 
obtaining  U.  Poor  results  are  obtained  because  of  random 
particle  noise  appearing  in  <uv> .  Values  of 
<u 2 > 1 7  2 / <U> s  =0 .16  and  <v 2 > 1 7 2/<U> s=0 . 23  were  obtained  at 
2400  time  steps  and  values  of  <u2 > 1 7 2/<Us >=0 . 1 7  and 
<v2> 1 7 2/<Us >=0 . 25  occurred  at  3600  time  steps,  with  25 
particles  per  cell.  With  9  particles  per  cell,  both 
turbulent  intensities  increased  in  time  from 
<u2> 1 7 2/<Us >=0 . 1 8 ,  <v2> 1 7 2/<Us >=0 . 1 7  at  600  time  steps  to 
0.25,  0.30  at  1800  time  steps.  Values  of  <u2 > 1 7 2/<Us >=0 . 09 
and  <v2> 1 7 2/<Us>=0 . 07  were  expected  from  the  discussion 
which  yielded  equation  (3.9A).  The  code  values  for 
<v 2  > 1 7  2 /<u  s  >  are  higher  than  <u2>172/<Us>  because  the 
Gaussian  averaging  of  fluid  quantities  was  only  done  in  the 
x-direction.  As  a  result,  noise  levels  were  higher  in  the 
y-direction  and  also  increased  more  rapidly  in  time. 
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However,  after  it  was  felt  that  the  channel  had  reached 
a  time  step  which  was  reasonably  steady  state  with  25 
particles  per  cell,  a  value  of  am=0.1  was  used  in  the 
y-momentum  equation  (equation  [A. 2])  for  the  next  600  time 
steps.  This  was  done  to  remove  the  large  random  noise 
levels  appearing  in  the  fluid  velocity  because  of  the  random 
density  fluctuations,  mentioned  in  section  4.5.  A  value  of 
am=0.0015  was  still  used  in  the  x-momentum  equation 
(equation  [ A . 1 ] )  and  a  value  of  at=0.5  was  still  used  in  the 
temperature  equation.  Within  a  few  hundred  time  steps, 
random  noise  levels  became  very  small  and  the  shear  stress 
profile,  mostly  determined  by  -p<uv>,  agreed  well  with  the 
known  profile,  as  is  seen  in  figure  4.7.  (Note:  The 
computational  cells  next  to  the  wall  extended  slightly 
beyond  the  viscous  sublayer  so  that  the  viscous  shear  stress 
pdU/9y  was  much  smaller  than  the  Reynolds  shear  stress 
-p<uv>  for  all  cells  except  those  next  to  the  wall,  as 
expected  [section  3.2].)  Values  of  <u2 > 1 7 2/<Us >=0 . 1 1  and 
<v2>  1  7  2/<Us>  =  0 . 07  were  observed  at  this  point,  which 
compares  favorably  with  the  expected  values  of  0.09  and 
0.07.  In  the  next  600  time  steps,  however,  using  such  a 
large  value  of  am=0.1  in  the  y-momentum  equation  interfered 
with  the  viscous  terms,  as  mentioned  earlier,  and  the  code 
progressed  to  an  incorrect  solution. 

Though  there  was  too  much  noise  in  <uv>  to  obtain  a 
proper  shear  stress  profile  with  values  of  am=0.0015  in  both 
the  x  and  y-momentum  equations  and  at=0.5;  it  is  clear  that 
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Figure  4.7  Comparison  of  code  shear  stress  profile 

(triangles)  to  known  profile  (straight  line). 
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the  actual  shear  caused  by  r  =  nd\J/dy  -  p< uv>  was  correct, 
since  a  correct  velocity  profile  is  observed  and  also  a 
correct  shear  stress  profile  is  seen  if  the  noise  is 
removed . 


CHAPTER  V 


Simulation  of  a  Magnetically  Stabilized  Gas  Discharge 

5.1  Analytic  considerations 

The  initial  purpose  of  this  project  was  to  develop  a 
viscous  fluid  particle  code  for  simulating  a  recently 
developed  method  of  stabilizing  a  gas  discharge,  which  was 
being  experimentally  tested  at  the  University  of  Alberta. 
Several  discharge  geometries  have  been  developed  at  the 
U  of  A,  however,  one  which  had  been  previously  modelled  by 
D.  Antoniuk  with  an  Alternating  Direction  Implicit  (ADI) 
fluid  code,  as  described  in  [4]  and  [16],  was  chosen  so  that 
the  fluid  particle  approach  could  be  compared  to  an  Eulerian 
fluid  approach  in  its  ability  to  model  a  gas  discharge. 

The  system  that  was  modelled  is  described  by 
D.  Antoniuk  in  [16],  and  in  more  detail  in  [4],  These 
descriptions  are  summarized  here  for  completeness.  The 
discharge  was  a  transverse  electrode  structure  which  was 
comprised  of  two  copper  plate  electrodes  with  a  d.c.  mag¬ 
netic  field  introduced  by  an  electromagnet  placed  at  the 
cathode,  as  depicted  in  figure  5.1.  The  magnetic  field  is 
shown  in  figures  5.2  and  5.3,  which  are  also  from  [4]  (note 
B©=0).  The  maximum  field  strength  in  figure  5.3  was  about 
0.2  Tesla  for  the  simulations  done  here.  The  gas  in  the 
discharge  was  a  mixture  of  helium,  nitrogen,  and  carbon 
dioxide.  As  described  in  [4]  and  [16],  the  total  magnetic 
field  (B)  of  the  discharge  can  be  approximated  very  well  as 
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Figure  5.1  Magnetically  stabilized  transverse  discharge 

geometry.  This  figure  is  taken  from  [4]  with 
the  author’s  permission. 
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Figure  5.2  Magnetic  flux  density  of  electromagnet  -  figure 


taken  from  [4]  with  author's  permission. 
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Figure  5.3  Magnetic  field 
electromagnet 


intensity  contours  of 
-  figure  taken  from  [4]  with 
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being  solely  due  to  the  electromagnet  since  virtually  no 
magnetic  induction  is  produced  by  the  current  flowing 
through  the  discharge.  The  current  density  flowing  between 
the  copper  plates  can  be  approximated  by  assuming  it  is 
axial,  i.e.  J©=Jr=0,  and  as  varying  only  in  the  radial 
direction.  As  in  [4],  a  current  density  which  was  constant 
in  time,  of  the  form: 

Jz(r)  =  J0exp[- (  r-r0  )  2/2b2  ]/b(27r )  1  7  2  (5.1) 

was  used;  where  r0  is  the  radial  distance  to  where  the  peak 
current  density  occurs;  b2=c2/8-^n2  where  c  is  the  distance 
between  the  points  where  Jz  is  one-half  its  maximum  value; 
and  J0  is  such  that  the  total  current  (It)  flowing  though 
the  discharge  (It  =  1 1 J z ( r ) rdrd# )  matches  the  current 
flowing  in  the  wires  to  the  electrodes.  This  profile  is 
shown  in  figure  5.4.  This  corresponds  to  an  annular  glow 
discharge,  which  was  observed  experimentally  by  Antoniuk. 

In  addition,  a  constant  value  of  electric  field,  E,  with 
only  an  axial  (z)  component  was  used. 

With  these  assumptions  made,  the  complete  MHD  equations 
simplify  considerably,  since  J,  B  and  E  are  no  longer 
variables.  Just  as  in  channel  flow,  the  fluid  quantities  v, 
p,  T,  and  p  are  the  unknowns,  and  with  only  slight 
modifications  to  equations  (2.2)  and  (2.6),  the  equations 
necessary  for  simulating  this  discharge  can  be  written.  It 
is  necessary  to  add  the  JXB  Lorentz  force  to  the  momentum 
equation  (2.2)  and  the  Joule  heating  term  J*E  to  the  tem¬ 
perature  equation  (2.6)  [17].  One  more  simplifying 
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Figure  5.4  Gaussian  shaped  current  density  distribution  - 

figure  taken  from  [4]  with  author's  permission 
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assumption  can  be  made  by  noting  that  the  Joule  heating  term 
will  dominate  the  viscous  dissipation  <p  from  equation  (2.6). 
This  is  because  the  ratio  of  viscous  dissipation  to  thermal 
heating,  defined  as  the  Eckert  number: 

E  =  V2/(Cp) (AT)  (5.2) 

(where  V  is  the  average  magnitude  of  the  total  fluid 
velocity;  Cp  is  the  specific  heat  at  constant  pressure;  and 
AT  is  the  characteristic  temperature  difference),  is  much 
less  than  one  for  the  gas  being  modelled  here.  The 
resulting  equations  for  the  magnetically  stabilized  gas 
discharge  are  written  in  cylindrical  coordinates  r  and  z  in 
Appendix  B.  Because  of  the  azimuthal  symmetry  involved,  6 
itself  was  not  a  variable,  however,  the  azimuthal  velocity 
ve  was  a  variable  i.e.  ve=v©(r,z),  resulting  in  a  so  called 
"2&1/2  dimensional"  model. 

5.2  Computational  considerations 


5.2.1  Boundary  and  initial  conditions 

The  computational  aspect  of  this  system  is  much  simpler 
than  the  channel  since  there  are  no  open  flow  boundaries. 

The  computational  window  was  chosen  as  illustrated  in  figure 
5.5.  Particles  were  reflected  from  the  anode  and  cathode 
walls.  In  addition,  at  r=0  the  particles  were  reflected 
since  this  is  a  line  of  symmetry.  The  radius,  R,  of  the 
outer  boundary  of  the  window  was  chosen  large  enough,  based 
on  the  results  of  Antoniuk,  so  that  for  r>R  the  state  of  the 
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Figure  5.5  Arrangement  of  ghost  point  boundaries  on  the 

finite  difference  mesh  -  figure  taken  from  [4] 
with  author’s  permission 
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gas  did  not  change  significantly  from  its  initial  condition. 
Thus,  near  the  outer  boundary  only  very  low  fluid  velocities 
were  expected  so  that  the  particles  were  reflected  from  the 
r=R  boundary;  this  kept  the  system  from  artificially  losing 
too  many  particles,  which  had  been  a  problem  with  the 
downstream  boundary  of  the  channel  before  pressure  clipping 
was  applied,  as  mentioned  in  Chapter  4.  The  fluid  boundary 
conditions  were  chosen  to  be  the  same  as  those  used  in  the 
ADI  simulation  of  [4]  and  [16],  in  order  that  the  results 
could  be  readily  compared.  On  the  r=0  boundary,  a  v=0  con¬ 
dition  was  applied  in  [4]  and  [16],  although  it  might  seem 
more  appropriate  to  use  a  Neumann  condition  for  the  axial 
velocity  v2  and  a  homogeneous  Dirichlet  condition  for  ve  and 
vr*  However,  to  be  consistent  with  [4]  and  [16],  v=0  was 
used  so  that  for  the  ghost  points  on  the  line  r=-A/2: 
v  i  i  =  — _v  i  *  ,  j .  Temperature  and  pressure  were  given  Neumann 
conditions  at  this  boundary.  The  ghost  points  on  z=-A/2 
satisfy  v  ,  ,  l  =-v  ,  ,  ■)  *  ,  ,  Ti(j  =  2T0  “  Ti(j  +  1,  Pi,j=Pi(j  +  i,  where 
j= 1 ,  so  that  on  z=0:  v=0,  T=T0 ,  3p/3z=0.  Similarly  the 
ghost  points  on  z=d+A/2  were  given  fluid  values  which 
yielded:  v=0,  T=T0 ,  3p/3z=0  for  z=d.  At  r=R  a  Neumann  con¬ 
dition  was  applied  for  pressure  and  temperature.  Again,  to 
be  consistent  with  [4],  the  condition  for  velocity  at  r=R 
was:  vr=0,  vz=0  and  3vG/3r=0. 

Initially,  the  gas  was  completely  at  rest  with  nine 
particles  per  cell,  and  all  particles  had  a  temperature  of 
T0  =  room  temperature.  For  simplicity,  the  phenomenological 
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parameters  p,  Cw  k.  ,  Cv  of  the  gas  were  assumed  constant, 
although  it  would  be  possible  to  use  some  known  functional 
form  (e.g.  p<xT1/2)  or  a  table  lookup  for  each  of  these,  how¬ 
ever  the  fluid  equations  (2.2)  and  (2.6)  would  have  to  be 
rederived  allowing  for  the  temporal  and  spatial  variation  of 
these  parameters.  The  simulated  gas  was  a  2/8/10  Torr  mix¬ 
ture  of  C02/N2/He.  The  initial  pressure  p0  was  thus  20 
Torr,  so  that  the  initial  density  p0  was  p0  =  p0/RT0.  (The 
R  used  here  is  the  ideal  gas  constant,  not  to  be  confused 
with  the  R  indicating  the  outer  radius  of  the  discharge.) 

The  individual  specific  heats  of  N2  and  C02  were  assumed 
constant  at  values  obtained  for  a  temperature  of  500°C  from 
[19].  Cv  for  He  was  taken  as  3R/2.  The  specific  heat  of 
the  mixture  was  then  obtained  by  a  mass  weighted 
(gravimetric)  average  of  the  individual  specific  heats.  The 
thermal  conductivity  of  the  mixture  was  estimated  as 
k=760X10'4  J/sm°C  for  T=500°C  based  on  tables  of  C02-N2  mix¬ 
tures  and  helium-other  gas  mixtures  from  [19].  (Note:  using 
a  gravimetric  weighting  for  the  thermal  conductivity  of  mix¬ 
tures  is  incorrect).  The  viscosity  was  chosen  as  p=6.6X10~‘ 
kg/ms  as  in  [4]  and  [16].  A  19X11  computational  grid  window 
was  used  with  grid  spacing  A=5.1X10'3m  so  that  R=19A,  d=11A. 
A  total  current  (It)  of  five  amps  was  used,  with  exactly  the 
same  current  density  distribution  as  shown  in  figure  5.4 
i.e.  ro=3.0X10"2  m,  c=2.5X10'2  m.  A  magnetic  field  which 
corresponded  to  a  coil  with  the  following  characteristics 
was  used:  inner  radius  of  windings  =  2.5  cm;  outer  radius  of 
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windings  =  5.5  cm;  length  of  coil  =  7.5  cm;  number  of 
turns  =  5600;  current  flowing  through  the  coil  =  2.0  A.  The 
magnetic  field  values  were  obtained  from  Antoniuk,  who  used 
the  computer  program  POTENT  [20]  to  solve  numerically  for 
the  magnetic  field  as  a  function  of  r  and  z.  An  electric 
field  of  Ez=-7130  V/m  was  used  based  on  an  experimental 
voltage  drop  of  400  V  across  the  discharge  [4].  The  total 
voltage  drop  between  the  electrodes  was  about  525  V  with 
It=5.0  A,  but  about  125  V  occurred  across  the  very  thin 
nonlinear  cathode  fall  region  right  next  to  the  cathode. 

This  region  was  ignored  for  the  purpose  of  this  simulation, 
since  it  does  not  obey  the  MHD  equations,  and  was  much 
smaller  than  a  grid  spacing.  At  was  chosen  as 
At=0.17A/CsO  =  2.0  microseconds,  where  C * 0= ( 7RT0 ) 1 '  2  .  The 
particles  were  given  a  value  of  a2=0.5A2,  and  the  Gaussian 
fluid  averaging  used  this  same  value  of  a2.  It  should  also 
be  mentioned  that  a  two  dimensional  Gaussian  averaging  was 
done  for  the  fluid  velocity  and  temperature  i.e.  the  sums  in 
equation  (4.1)  were  extended  to  be  over  all  grid  points,  not 
just  over  the  cells  in  the  jth  row.  Note  that  with 
a2=0.5A2,  the  averaging  for  cell  i,j  is  essentially  only 
over  cells  i , j ± 1  and  i±1,j,  since  f  from  equation  (2.8)  is 
0.37  for  these  cells;  while  for  i±1,j±1:  f  =  0 .  14,  and  for 
i±2,j  or  i ,  j ± 2 :  f  =  0.02.  A  lower  value  of  a2/A2  was  used 
since  it  was  found  that  noise  levels  were  low  enough  with 
this  value.  (Note:  without  the  Gaussian  averaging,  random 
noise  levels  about  doubled) 
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5.2.2  Finite  differencing  considerations 

All  terms  in  the  MHD  equivalent  to  equations  (2.2)  and 
(2.6)  were  differenced  either  flux  or  force  conservatively 
except  the  v«(V«7r)  term  of  the  temperature  equation.  As 
with  the  channel,  the  additional  computational  effort 
involved  in  evaluating  flux  and  force  differencing  for  this 
term  was  not  deemed  worthwhile,  since  its  effect  (viscous 
work)  was  not  of  primary  interest.  Note  that  in  the  cells 
adjacent  to  the  r=0  boundary,  force  conservative 
differencing  cannot  be  used  on  terms  of  the  form  (9F/3r)/r 
since  a  force  conservative  differencing  would  yield: 
[(9F/9r)/r], . j  =  [(Fi,j-Fi_1(j)/ri_1/2,j 

+  (F  i  (  j  +  i“Fj  (  j)/ri  +  1/  2,  j ]/2A 
but  ri.1/2lj=0  for  these  cells,  so  that  1/ri_1/2(j  is 
undefined.  For  such  terms  where  F  was  a  single  variable, 

(F  j  +  i  (  j -F  j  .  ■,  (  j  )  /r  j  ;  j  A  was  used.  If  F  =  rG,  where  G=G(v,p,T)  , 
then  the  term  was  differentiated  i.e. 

(9F/9r)/r  =  G/r  +  9G/9r  where  3G/3r  was  differenced  flux 
conservatively,  and  G/r  =  G| ,  j/ri , j. 

5.3  Results  and  discussion 

A  run  of  2000  time  steps  (4  ms)  was  done  using  the 
boundary  and  initial  conditions  mentioned  earlier.  A 
contour  plot  of  vG  is  shown  in  figure  5.6.  Previous  time 
steps  have  very  similar  contour  plots  for  v©  except  that  the 
value  of  v©  increases  in  time,  so  that  they  are  essentially 
scaled  versions  of  that  shown  in  figure  5.6. 
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Figure  5.6  Contour  plot  of  ve  in  m/s  from  particle  code  at 

t=4  ms  (2000  time  steps). 
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Unfortunately,  because  we  have  assumed  6  symmetry,  any 
Rayleigh-Taylor  instabilities  [17]  in  6  which  might  arise  in 
the  actual  discharge  cannot  occur  in  the  simulation.  What 
is  meant  by  a  Rayleigh-Taylor  instability  is  shown  in  figure 
5.7.  The  initial  position  of  a  section  of  plasma  is  indi¬ 
cated  by  the  circle.  The  circumstances  under  which  this 
instability  will  occur  are  those  where  a  heavier  fluid  is 
" supported"  by  a  lighter  fluid.  In  this  case,  a  hot,  lower 
density  region  is  assumed  to  be  in  equilibrium  with  a 
denser,  lower  temperature  region,  while  a  centrifugal  force 
pushes  the  heavier  fluid  onto  the  lighter  one.  The 
instability  causes  portions  of  the  discharge  to  flow 
radially  outward  (indicated  by  the  arrows  pointing  radially 
outward)  while  neighboring  sections  flow  radially  inward. 
This  tends  to  mix  the  plasma  so  that  any  variations  in  the 
state  of  the  plasma  are  diffused.  The  wavy  line  indicates 
the  position  of  the  same  plasma  mass  shortly  after  the  onset 
of  the  instability.  Note  also  that  no  6  mixing  or  6 
turbulence  is  modelled  in  a  2&1/2  dimensional  model,  so  that 
the  computational  system  is  artificially  too  "calm". 

In  our  simulation,  the  Joule  heating  term  J2EZ  produces 
a  temperature  profile  which  is  reasonably  similar  in  shape 
to  that  of  J2  (since  E2  is  a  constant).  Such  a  temperature 
profile  coupled  with  the  centrifugal  force  resulting  from  v© 
would  certainly  cause  Rayleigh-Talor  instabilities  to  form; 
but,  with  just  the  two  spatial  dimensions  r  and  z,  it  is 
impossible  to  simulate  these  instabilities,  so  that  the 
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Figure  5.7  A  Rayleigh-Taylor  instability,  as  shown,  would 

normally  result  in  the  gas  discharge  geometry 
where  a  less  dense  fluid  contains  a  more  dense 
fluid.  See  p.  95  for  an  explanation. 
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Joule  heating  term  would  heat  a  ring  of  plasma  to  very  high 
temperatures,  while  the  r=R  wall  is  held  at  room  tempera¬ 
ture.  Clearly,  this  is  not  a  physical  situation.  Thus,  to 
obtain  a  useful  result  from  our  2&1/2  dimensional 
simulation,  it  is  necessary  to  stop  the  code  after  a  number 
of  time  steps  and  use  the  results  of  this  time  step  as  an 
approximation  to  the  steady  state  results. 

It  was  found  that  if  a  simulation  was  stopped  after 
tmx  =  4  ms,  where  tmx  is  the  time  required  for  the  average  gas 
temperature  to  increase  from  its  initial  value  of  room  tem¬ 
perature  to  a  value  in  the  range  of  experimental  values  from 
[4],  then  vG  (which  at  t=0  was  zero)  was  also  in  the  range 
of  experimental  data.  The  critical  variable  for  which  good 
experimental  data  exists  [4]  is  the  peak  value  of  ve  (vGmx). 
It  was  also  found  that  if  temperature  was  not  used  as  a 
variable,  then  the  amount  of  time  required  for  vGmx  to  reach 
a  steady  state  value  was  also  about  4  ms. 

The  particle  code  results  correspond  very  well  with 
those  obtained  by  D.  Antoniuk  both  experimentally  and 
computationally  in  [4].  Antoniuk' s  code  was  an  isothermal 
ADI  (Alternating  Direction  Implicit)  Eulerian  fluid  code 
which  assumed  a  constant  temperature  of  425  K.  His 
computational  results  indicate  a  steady  state  value  of 
v©mx=35  m/s  occurring  after  4  ms.  A  contour  plot  of  the 
steady  state  (t=4  ms)  sheared  rotational  velocity  distri¬ 
bution  produced  by  his  ADI  program  is  shown  in  figure  5.8. 
Experimentally,  Antoniuk  observed  a  steady  state  value  of 
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Figure 


5.8  Contour  plot  of  v©  at 
[16]  using  ADI  fluid 


t=4  ms  obtained  by  Antoniuk 
code . 
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vemx=38.5  m/s  [4,  p.  168]  with  a  total  discharge  current  of 
It=5.0  A  and  a  magnetic  field  associated  with  the  magnetic 
coil  current  Im=2.0  A.  (Note  that  in  [4],  two  equivalent 
coils  were  used;  one  required  2.0  A,  as  mentioned  here,  to 
give  the  same  field  used  in  the  particle  code,  while  the 
other  required  4.0  A  to  give  the  same  field.  The  data  on 
[4,  p.  168]  for  vemx  uses  the  4.0  A  coil.  All  Im  values 
listed  in  the  work  herein  are  for  the  2.0  A  coil.) 

Antoniuk  recorded  maximum  values  of  vr  in  the  range  of 
10  to  20  m/s  [4,  p.  174-177]  for  a  1/7/12  C02/N2/He  mixture. 
A  typical  velocity  flow  plot  of  fluid  velocity  Vj  }  (from 
equation  [4.1])  is  shown  in  figure  5.9,  where  the  maximum 
velocity  is  36.4  m/s.  Antoniuk  also  found  that  a  circular 
eddy  developed  in  his  simulation.  This  eddy  covered  the 
entire  computational  window  in  the  r-z  plane  of  his 
isothermal  ADI  simulation.  Although  secondary  flow 
developed  with  the  fluid  particle  simulation,  no  such  large 
scale  vortex  developed.  Isothermal  runs  were  done  to  deter¬ 
mine  if  the  temperature  profile  stopped  the  vortex  from 
forming,  but  no  such  large  scale  vortex  developed  in  either 
isothermal  or  nonisothermal  runs.  It  is  difficult  to  say 
whether  such  a  vortex  was  observed  experimentally.  Since 
Antoniuk' s  ADI  simulation  had  difficulty  dealing  with 
reversals  in  flow  velocities  occurring  between  cells,  it  is 
possible  that  the  vortex  structure  was  a  computationally 
preferred  flow  mode  of  the  ADI  simulation. 
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Figure  5.9  Typical  particle  code  velocity  flow  plot.  The 
arrows  length  are  scaled  with  the  magnitude  of 
velocity.  The  maximum  veloctiy  is  36.4  m/s. 
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It  is  difficult  to  compare  any  other  code  results  to 
Antoniuk's  results  because  most  of  Antoniuk's  data  was 
obtained  for  considerably  different  discharge  conditions. 

He  obtains  a  temperature  profile  [4,  p.  217]  for  a  discharge 
current  of  I  ,*4.0  A,  but  using  a  cathode  which  was  composed 
of  many  metal  pins,  not  a  flat  plate.  The  multipin  cathode 
causes  the  temperature  profile  to  be  less  peaked  than  for 
the  flat  plate  cathode.  In  addition,  a  discharge  current  of 
I ,=4.0  A  produces  a  discharge  with  a  current  density  that 
has  its  peak  value  in  a  different  place  than  that  in  figure 
5.4.  As  well,  there  will  be  less  heating  with  only  4.0  A 
flowing  through  the  discharge.  Even  so,  one  can  estimate 
that  the  temperatures  obtained  from  the  code  at  t=4  ms,  as 
in  figure  5.10,  are  quite  reasonable,  as  mentioned  earlier. 

It  should  be  noted  that  the  value  of  mult i streaming 
coefficient  am  used  for  this  geometry  was  very  much  larger 
than  that  used  in  the  channel.  The  value  of  am=0.1  used 
here  can  produce  artificial  drag  forces  which  are  up  to 
several  thousand  times  larger  than  the  viscous  forces  caused 
by  m.  (An  analysis  similar  to  that  which  yielded  equation 
[4.9]  yields  a  value  of  am=0.0001  here.)  This  is  in 
contrast  to  the  channel  where  it  was  necessary  to  set  am  so 
that  artificial  drag  forces  and  viscous  forces  were  of  the 
same  order.  As  mentioned  in  chapter  4,  it  is  not  necessary 
to  restrict  am  in  this  way  since  there  is  no  boundary  layer 
flow  in  the  r-z  plane  and  the  flow  in  the  6  direction  is 
artificially  laminar  because  the  code  is  two  dimensional. 
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Figure  5.10  Contour  plot  of  gas  temperature  T  in  degrees 

Kelvin  from  particle  code  at  t=4  ms. 
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Thus,  there  is  no  need  to  have  intracell  and  intercell 
viscous  effects  compatible  so  that  there  is  no  difficulty  in 
using  a  value  similar  to  the  value  of  at=am=0.5,  used  by 
Leboeuf  et  al.  in  [2]  for  an  inviscid  fluid. 

It  should  also  be  mentioned  that  if  a  conventional  Vp 
finite  differencing  was  used,  values  of  the  turbulent 
kinetic  energy  <q2>1/2  (as  in  equation  [3.10])  increased  by 
a  factor  of  ten  or  more  from  the  modified  Vp  case  to  levels 
that  produced  very  unrealistic  results.  This  increase  in 
<q2>1  2  was  due  to  large  scale  random  motions,  as  mentioned 
in  section  2.2.2.  1,  and  was  only  kept  stable  by  the  Gaussian 
averaging  combined  with  the  mul t i st reami ng  term  a.  Using 
the  modified  Vp  differencing  and  values  of  am=0.1,  at=0.5, 
random  particle  noise  levels  were  small. 


CHAPTER  VI 


Concluding  Remarks 


6 . 1  Summary 

An  approach  for  applying  a  NGP  (nearest  grid  point) 
finite  particle  method  to  viscous,  compressible  fluids  has 
been  described.  The  fluid  is  modelled  very  simply  by  a  set 
of  discrete  finite  particles.  The  particles  carry  a  value 
of  velocity  and  temperature,  which  give  rise  to  fluid 
velocities  and  temperatures  on  an  Eulerian  grid  by  averaging 
the  particle  velocities  and  temperatures  in  each  cell  on  the 
grid.  The  particles  also  have  a  Gaussian  mass  density 
shape,  allowing  the  fluid  density  to  be  calculated  by 
summing  the  particle  densities  using  a  NGP  approach.  The 
instantaneous  fluid  equations  were  written  in  a  Lagrangian 
conservative  form  so  that  the  particle  temperatures  and 
velocities  were  pushed  conservatively  in  time  using  the 
Lagrangian  derivatives  of  velocity  and  temperature,  which 
were  known  by  finite  differencing  the  fluid  equations.  The 
Vp  term  in  the  velocity  equation  was  dealt  with  in  a  special 
manner  which  removed  the  previous  difficulty  of  excessively 
high  random  particle  noise  levels  appearing  with  the  NGP 
approach.  In  the  channel  flow  simulation,  a  small  value  for 
the  multistreaming  coefficient  am  (appearing  in  the  momentum 
equation)  was  necessary  otherwise  an  incorrect  velocity 
profile  developed.  Analytic  and  experimental  results  for 
channel  flow  were  illustrated.  A  particle  model  was 
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developed  for  turbulent  channel  flow,  with  the  use  of  buffer 
zones  at  the  open  ends  and  pressure  clipping  at  the 
downstream  end  to  maintain  a  steady  flow  rate.  The 
constraints  on  various  parameters  necessary  for  a  conver¬ 
gent,  stable  solution  were  developed.  To  maintain 
stability,  it  was  found  that  the  temperature  equation  was 
required  for  channel  flow.  Excellent  agreement  with  experi¬ 
mental  results  for  the  mean  value  of  pressure  and  along 
stream  velocity  for  turbulent  channel  flow  (Reynolds 
number=5333)  were  obtained.  Profiles  of  the  turbulent  cor¬ 
relation  <uv>  agreed  well  with  analytic  results  if  the 
random  noise  was  reduced  to  low  levels,  otherwise  <uv>  was 
covered  by  noise  which  was  implicit  to  this  average.  The 
noise  however  did  not  cause  turbulent  transport  of  momentum, 
as  seen  by  the  correct  velocity  profile  and  correct  shear 
stress  profile. 

A  magnetically  stabilized  gas  discharge  was  also 
modelled.  Good  results  were  obtained  in  comparison  with  the 
experimental  data  of  Antoniuk  [4]  and  with  the  computational 
results  of  Antoniuk' s  isothermal  ADI  conventional  Eulerian 
fluid  code. 

6.2  Some  further  remarks 

The  NGP  particle  approach  to  modelling  a  fluid  is  a 
very  straightforward  method  that  is  easily  understood  and 
which  is  also  very  physical  in  how  it  deals  with  the  fluid. 
As  a  result,  difficulties  with  the  method  can  usually  be 
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understood  on  a  very  practical  level  without  requiring  much 
detailed  analysis.  Because  each  cell  has  many  particles, 
the  resolution  of  a  flow  is  also  greater  than  that  of  con¬ 
ventional  fluid  methods. 

A  typical  run  for  the  channel  required  about  75  (89) 
minutes  of  computing  time  per  600  time  steps,  for  9  (25) 
particles  per  cell,  but  85%  (70%)  of  this  computing  time  was 
taken  up  by  the  subroutines  used  for  calculating  the 
Gaussian  densities  and  Gaussian  averaging.  Because  these 
were  done  on  the  array  processor  (FPS  190L),  which  cost  1/20 
that  of  the  mainframe  computer  (Amdahl  5860  which  is  equiv¬ 
alent  to  an  IBM  370),  the  cost  of  modelling  the  gas 
discharge  was  only  one  and  one  half  to  two  times  that  of 
Antoniuk's  ADI  simulation.  With  the  gas  discharge,  a 
typical  run  required  about  seven  minutes  mainframe  CPU  time 
and  18  minutes  of  array  processor  time. 

6.3  Future  applications  and  considerations 

Clearly,  there  is  much  room  for  future  applications  of 
the  fluid  particle  method.  The  work  herein  has  developed 
necessary  considerations  when  dealing  with  turbulent  viscous 
compressible  fluids  in  both  open  and  closed  systems  using  a 
fluid  particle  approach.  It  would  be  most  interesting  to 
determine  how  closely  more  of  the  details  of  actual 
turbulence  are  modelled  by  a  fluid  particle  code.  It  would 
of  course  be  necessary  to  remove  the  random  noise  from  the 
turbulence.  If  enough  particles  are  used  per  cell  in  a  NGP 
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approach,  the  random  particle  noise  appearing  for  small  am 
should  be  reduced  to  low  levels,  as  indicated  by  equations 
[4.15]  and  [4.16].  Because  the  Gaussian  subroutines  take 
more  than  70%  of  the  computing  time  and  because  they  depend 
on  the  number  of  cells,  and  not  the  number  of  particles  per 
cell  (n0),  the  total  computing  time  increases  by  only  a  few 
percent  by  doubling  the  number  of  particles  per  cell. 
Unfortunately,  at  the  U  of  A,  this  increase  was  in  mainframe 
CPU  time  and  not  AP  time,  so  that  costs  increased  by  a 
factor  of  about  1.5  for  n0  doubling.  On  a  vector  processing 
supercomputer  (e.g.  a  Cyber  205  or  Cray  1)  one  would  likely 
put  the  whole  program  on  the  supercomputer  so  that  runs  with 
large  n0  would  require  only  a  small  percentage  more  time. 

It  is  felt  that  noise  levels  in  the  channel  would  likely  be 
reduced  to  small  levels  if  at  least  81  particles  per  cell 
were  used. 

Other  methods  of  grid  particle  assignment  might  be 
attempted  to  reduce  noise  levels,  but  these  would  likely 
require  many  special  considerations  unique  to  each  method. 
The  noise  levels  in  the  area  weighting  method  listed  in  [7] 
would  have  to  be  reduced  before  being  usable  in  turbulent 
flow,  while  the  subtracted  dipole  approach  (SUD)  described 
in  [7]  would  require  special  modification  for  nonperiodic 
flow  because  it  uses  fast  Fourier  transforms  (FFT)  for  the 
spatial  derivatives,  so  that  leakage  and  possibly  aliasing 
errors  would  appear. 


108 


Areas  for  further  work  include:  the  simulation  of  a 
viscous  flow  past  an  airfoil;  extension  of  the  fluid 
particle  method  to  incompressible  fluids;  the  development  of 
a  NGP  approach  with  variable  grid  spacings;  the  use  of 
coordinate  transformation  techniques  to  simulate  irregular 
geometries  or  to  transform  a  variably  spaced  grid  into  a 
regularly  spaced  one  and  application  of  a  fluid  particle 
approach  to  this  set  of  equations. 

In  conclusion,  as  Hockney  and  Eastwood  comment  [21, 
p.  347],  "the  principle  of  moving  the  fluid  as  particles  in 
order  to  ensure  positive  density  and  mass  conservation,  and 
calculating  the  force  on  the  fluid  particles  from  field 
quantities  derived  on  a  fixed  mesh,  is  sound  and  likely  to 
have  wide  application  in  fluid  dynamics." 
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Appendix  A 

The  full  compressible  fluid  equations  (2.2)  and  (2.5)  in  two 
dimensional  cartesian  coordinates  x  and  y  are: 
pdvx/dt  =  -9p/9x  +  p[32vx/9x2  +  dJvx/9y2] 

+  p[92vx/9x2  +  92vy/9x9y]/3  (A.1) 

pdvy/dt  =  -9p/9y  +  p[92vy/9x2  +  92vy/9y2] 

+  p[92vx/9x9y  +  92vy/9y2]/3  (A. 2) 

pCv(dT/dt)  =  k[92T/9x2  +  92T/9y2]  -  9(pvx)/9x  -  9(pvy)/9y  + 
^{vx ( 9 2 vx/9x 2+9 2vx/9y 2 )  +  v y ( 9 2 v y /9x 2 + 9 2 v y /9y 2 )  + 

[  vx  (  9  2  vx/9x  2 +  9  2  vy /9x9y  )  +  v y  (  9  2  v y /9y  2  -t-9  2  v x /9x 9y  )  ] /3  + 

4 [ ( 9vx/9x+9vy /9y ) 2  -  ( 9vy /9y ) ( 9vx/9x ) ]/3  + 

( 9vy/9x+9vx/9y ) 2 }  -  pvxdvx/dt  -  pvydvv/dt  (A. 3) 

The  viscous  term  in  the  temperature  equation  (A. 3)  can  be 
written  in  a  flux  conservative  form  (see  Chapter  2),  however 
the  rewards  of  such  a  form  were  not  believed  to  be  worth  the 
excess  complexity  added  to  the  finite  differencing  since  the 
viscous  work  and  dissipation  were  not  of  major  concern,  as 
mentioned  in  Chapter  2. 
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Appendix  B 


The  full  compressible  fluid  equations  (2.2)  and  (2.6) 

written  for  two  cylindrical  coordinates  r  and  z  and  three 

velocity  components  v  r  v©,  and  vz  are: 

pdvr /dt  =  -3p/3r  +  p[  3j_r  3vr  /3r  )  ]  +  p[32vr/3z 2  -  vr  /r  2  3 

r  dr 

♦  p3_[_  <rvr  )/dr]/3  +  p ( 3 2 v2 /3z3r ) /3 
9r  r 

+  ( J©BZ-JZB© )  +  pve  2/ r  (B. 1 ) 

pdv©/dt  =  p[  3J_r  3v©/3r  )  ]  +  p[32v©/3z2  -  v8/rJ] 
r  9r 

+  (JzBr~JrB2)  -  pvevr/r  (  B  .  2  ) 

pdv2/dt  =  -9p/9z  +  p [  3J_r 9v z /9 r  )  ]  +  ( 4p3  2  v z /9z  2  ) /3 

r  9r 

+  J_p3  2  (  rvr  ) /9r  3z  +  (JrB©-J©Br)  (B.3) 

3r 

pCv(dT/dt)  =  -pvrdvr/dt  -  pv©dv©/dt  -  pvzdvz/dt 

•+  k  9  (  r  9T/9r  )  +  k32T/3z2 
r  9r 

-  {  [  9j_rpvr  )  ]/r  +  9(pvz)/9z] 

3  r 

+  p[vr(V2v)r  +  v© ( V2 v ) ©  +  v  z ( V2 v ) 2 
+  (vr32vz/3z3r  +  vz92vr/9z9r  -  vr 3 2 vr /3z 2 ) /3 
+  (vz 3vr/3z )/r ] 

+  JrEr  +  J©E©  +  JZEZ  (B.4) 

where  V2v  =  (V2v)fer  +  (V2v)©e©  +  (V2v)zez  and  the  e's 
are  unit  vectors. 


113 


